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EXECUTIVE  SUMMARY 


Communication  in  the  underwater  battle  space  is  essential  to  Netcentric  Warfare 
concepts  in  the  21st  century.  Acoustic  communication  holds  promising  capabilities  for 
connecting  the  undersea  environment.  The  underwater  acoustic  sound  channel  however, 
presents  challenging  obstacles  that  must  be  overcome  to  allow  satisfactory  acoustic 
communication.  Possibly  the  most  significant  are  the  frequency  selective  fading  and 
multipath  natures  the  underwater  channel  presents. 

Orthogonal  Frequency  Division  Multiplexing  (OFDM)  is  extremely  robust  to 
multipath  frequency  selective  fading  channels  and  therefore  is  rapidly  developing  interest 
in  application  to  underwater  acoustic  communication.  The  primary  thrust  of  this  work  is 
to  develop  a  computer-based  simulation  of  an  OFDM  based  communication  system  for 
the  undersea  environment.  The  result  will  then  provide  a  foundation  upon  which  the  key 
parameters  and  capabilities  of  underwater  acoustic  communication  using  OFDM  can  be 
investigated.  The  immediate  benefit  of  the  work  is  the  establishment  of  the  achievable 
communication  system  parameters,  common  to  all  communication  systems,  which  are  bit 
rate,  bit  error  rate,  bandwidth  and  the  required  signal  to  noise  ratios. 

The  simulation  is  performed  in  the  Matlab  environment  where  all  major  hardware 
in  the  transmitter  and  receiver  are  simulated.  The  underwater  acoustic  sound  channel  is 
modeled  using  the  physics  based  parabolic  equation  approximation.  The  model  is  the 
Monterey-Miami  Parabolic  Equation  (MMPE)  which  employs  the  Split  Step  Fourier 
method.  Use  of  a  physics  based  model  vice  the  typical  stochastic  model  for  the  acoustic 
channel  ensures  a  better  realization  of  the  phenomenon  that  is  experienced  in  the 
undersea  environment. 

The  simulation  allows  for  any  type  of  data  to  be  transmitted  from  analog  voice 
signals  to  binary  data.  Significant  flexibility  is  coded  into  the  simulation  to  allow  for  a 
wide  range  of  parameter  analysis  and  sensitivity  testing.  Reed  Solomon  coding  and  block 
interleaving  are  applied  to  achieve  several  decibels  of  forward  error  correction  coding 
gain.  Quadrature  amplitude  modulation  (QAM)  is  applied  with  constellation  sizes  of  8  to 
64  points.  Double  sideband  suppressed  carrier  modulation  is  used  to  modulate  the 
baseband  signal  to  the  desired  transmission  band. 

Bit  rates  achieved  range  from  4  to  10  kilobits  per  second  (kbps)  utilizing  from  2 
to  6  kilohertz  (kHz)  of  bandwidth  depending  on  the  size  of  the  QAM  constellation 
applied.  The  transmitted  signals  occupy  frequency  bands  between  4  and  14  kHz. 

The  results  are  validated  through  a  theoretical  development  of  achievable 
performance.  In  addition  validation  is  made  through  comparison  to  in  water  testing  of 
OFDM. 

The  results  substantiate  the  feasibility  of  OFDM  as  an  attractive  candidate  for 
underwater  acoustic  communication.  The  product  of  the  work  provides  validated  system 
performance  expectations  as  well  as  a  foundation  upon  which  future  work  in  this  area  can 
be  perfonned. 
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I.  INTRODUCTION 


A.  PURPOSE 

Netcentric  Warfare  has  moved  to  the  forefront  of  the  Navy’s  priorities  for  the  21st 
century.  Essential  to  the  concept  of  Netcentric  Warfare  is  the  ability  to  communicate 
quickly,  reliably,  efficiently  and  in  many  cases  covertly  in  all  regions  of  the  battle  space. 
More  specifically,  key  characteristics  of  the  communication  capability  required  are  high 
data  rates,  low  error  rates,  environmental  robustness  and  minimally  overt  acoustic 
signatures.  The  undersea  environment  is  arguably  the  most  challenging  region  of  the 
battle  space  in  which  to  communicate  effectively.  Current  and  future  goals  for 
underwater  wireless  communication  involve  data  flow  between  submarines,  surface  units, 
unmanned  undersea  vehicles  (UUVs)  and  ocean  surface  and  ocean  bottom 
communication  buoys  in  any  combination  via  the  underwater  sound  channel.  Research 
and  development  in  the  fundamental  communication  methods  and  protocols  that  support 
Netcentric  Warfare  concepts  will  enhance  the  Navy’s  ability  to  fight  successfully  in  all 
areas  of  the  battle  space. 

Orthogonal  Frequency  Division  Multiplexing  (OFDM)  possesses  characteristics 
that  have  the  potential  to  overcome  the  adverse  effects  of  frequency  selective  multipath 
fading  present  in  the  underwater  sound  channel.  Establishing  groundwork  in  this  area  to 
address  the  feasibility  of  OFDM  for  underwater  communication  is  essential  to  developing 
an  underwater  communication  system  that  is  beneficial  to  the  Navy. 


B.  GOALS 

The  primary  goal  is  to  develop  a  computer  simulation  that  will  allow  for  research 
into  the  key  parameters,  capabilities,  and  limitations  of  OFDM  in  underwater  acoustic 
communication.  The  simulation  provides  a  ‘test  bed’  for  future  study  of  other  modulation 
methods,  forward  error  correction  coding,  quantization,  signal  detection,  decoding  etc.  In 
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developing  and  validating  the  simulation  the  feasibility  of  OFDM  systems  in  underwater 
communication  is  verified.  There  is  one  known  experiment  with  OFDM  in  the 
underwater  acoustic  sound  channel  that  is  used  as  a  standard  against  which  the  results  of 
the  simulation  are  validated.  The  work  was  done  by  Coatelan  and  Glavieux.  [Ref.  1] 

While  there  has  been  considerable  work  using  OFDM  for  communication  all  but 
the  work  of  Coatelan  and  Glavieux  [Ref.  1]  and  that  of  Kim  and  Lu  [Ref.  2]  has  been 
outside  the  area  of  underwater  acoustic  communication.  Kim  and  Lu  [Ref.  2]  present  a 
new  broadband  underwater  channel  model  that  has  the  capability  to  model  the  Doppler 
time  scaling  effect  as  well  as  multi-path  effects.  The  propagation  path  delay  is  modeled 
with  two  components  that  generate  the  simulated  Doppler  shift  determined  by  the  relative 
motion  between  the  source  and  receiver.  The  simulation  work  by  Kirn  and  Lu  [Ref.  2] 
represents  the  only  previous  work  in  simulating  an  underwater  acoustic  communication 
system  employing  OFDM  as  the  backbone.  A  physics  based  model  is  used  to  provide  a 
realistic  model  of  underwater  acoustic  sound  channel. 

C.  METHODOLOGY 

The  simulation  is  run  on  a  PC  using  Matlab  code.  The  underwater  sound  channel 
is  modeled  using  a  physics  based,  parabolic  equation  (PE)  model,  namely  the  Monterey 
Miami  Parabolic  Equation  model. 

The  simulation  is  used  to  conduct  experimentation  to  examine  key  parameters 
such  as  achievable  bit  rate  and  bit  error  rate,  required  SNR,  transmission  range, 
robustness  to  multipath  propagation  and  frequency  selective  fading,  and  bandwidth 
requirements.  The  strengths  and  weaknesses  of  OFDM  applied  to  underwater  acoustic 
communications  will  be  identified  to  assess  its  feasibility  and  merit  for  future  work. 


D.  BENEFITS 

Simulation  is  rapidly  becoming  a  popular  method  for  experimentation  with 
technology-based  systems.  One  of  the  biggest  advantages  of  simulation  is  the  ability  to 
perform  countless  numbers  of  experiments  under  controlled  conditions  at  a  minimal  cost. 
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Some  of  the  biggest  obstacles  to  simulation  are  the  decision  of  what  level  of  simulation  is 
required  to  accurately  capture  the  real  world  and  validation  of  the  simulation.  With  all 
simulations  there  must  be  an  initial  model.  This  model  suffers  from  some  simplifications 
that  cause  departure  from  reality.  Later  improvements  will  help  the  model  to  converge  on 
the  characteristics  of  the  real  world  system  that  is  being  modeled.  With  that  in  mind  many 
design  decisions  have  been  made  in  developing  the  simulation.  The  emphasis  is  on 
modeling  the  aspects  of  the  system  which  have  the  greatest  impact  on  the  key  parameters 
of  interest. 

Finally,  the  benefit  of  this  thesis  is  the  ability  to  experiment  with  OFDM  based 
communication  systems  such  that  the  factors  influencing  the  key  parameters  in 
communication  can  be  understood.  Understanding  these  factors  will  provide  the 
knowledge  necessary  to  build  actual  systems  for  underwater  acoustic  communication. 
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II.  OFDM  THEORY 


A.  OFDM  BASICS 

Modulating  an  orthogonal  set  of  subcarriers  with  Q-ary  complex  QAM  symbols 
generates  an  OFDM  signal.  The  basic  OFDM  system  is  shown  in  Fig.  2.1  below.  The 
input  bit  stream  is  buffered  and  converted  to  parallel  complex  QAM  symbols  in  the  QAM 
Modulator. 


Input  Bit 
Stream,  R 

(bps) 


Output  Bit 
Stream 


Figure  2.1.  OFDM  Basic  Block  Diagram. 


The  ith  subchannel,  denoted  as  Xy<,  where  i  is  the  subchannel  index  and  k  is  the 


OFDM  symbol  index,  can  be  independently  modulated  as  in  adaptive  modulation 
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schemes  or  all  subcarriers  may  be  modulated  in  the  same  manner.  The  parallel  QAM 
symbols  are  then  applied  to  the  IFFT,  which  generates  N  parallel  real  samples.  The  real 
samples  are  then  converted  to  serial  format.  In  this  simplified  block  diagram  the  receiver 
reverses  the  operations  performed  by  the  transmitter. 

The  QAM  modulator  buffers  the  input  bits  into  blocks  of  size  b,  where  b  =  RT,  T 
is  the  symbol  duration  in  seconds  and  R  is  the  overall  bit  rate  in  bits  per  second.  The 
symbol  rate  then  is  Rs  =  1/T  symbols  per  second.  The  number  of  bits  per  subchannel 
determines  the  size  of  the  QAM  constellation,  Q.  For  example  with  q  =  4  bits  per 
subchannel,  Q  =  24,  which  is  16  QAM.  In  OFDM  the  frequency  spectrum  of  the  signal  is 
created  in  the  QAM  Modulator  and  then  the  time  domain  signal  results  from  passing  the 
output  through  the  IFFT  operation.  The  complex  QAM  symbols  derived  from  the  binary 
code  words  are  sequentially  placed  in  a  subcarrier  band  and  OFDM  symbol.  Figure  2.2 
illustrates  the  organization  of  four  OFDM  blocks  with  eight  subcarriers. 
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Figure  2.2.  Subcarrier  and  OFDM  Symbol  Arrangement. 


The  IFFT  provides  an  orthogonal  transformation  and  maintains  the  energy 
contained  in  the  input  symbol  [Ref.  3]  Written  as  an  equation,  this  is 


j  N  2 

N 

=2X* 

n- 1 

Due  to  the  orthogonality  of  the  IFFT  the  subcarriers  each  have  an  integer  number 
of  cycles  within  the  FFT  interval.  The  number  of  cycles  in  the  FFT  interval  for  a  given 
subcarrier  differs  by  one  from  that  of  adjacent  subcarriers.  Note  that  the  terms  subchannel 
and  subcarrier  are  used  interchangeably  in  OFDM  tenninology.  Figure  2.3  illustrates  the 
orthogonality  of  3  subcarriers  in  an  OFDM  symbol.  In  this  example  each  of  the 
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subcarriers  has  the  same  magnitude  and  phase  but  this  is  not  in  general  the  case.  As  can 
be  seen,  each  subcarrier  has  an  integer  number  of  cycles  within  the  FFT  interval. 


1.  Subcarrier  Orthogonality 

The  real  bandpass  OFDM  symbol  can  be  expressed  as  follows: 


2  /  I]  /  2 

40  =  Re(X  d  -j  exp(y'2;r(/c  - — )(t-t0))},  t0<t<t0+T 

,=_N  /2  1 

‘  2 

x(t )  =  0 ,t  <t0  t  >  t0  +  T 


(2.2) 


where  T  is  the  symbol  period,/^  is  the  carrier  frequency,  N  is  the  number  of  subcarriers, 
di  are  the  complex  QAM  symbols  and  to  is  the  symbol  starting  time  [Ref.  4].  From  Eq. 
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2.2  it  can  be  seen  that  the  OFDM  symbol  consists  of  non-zero  subcarriers  over  the 
symbol  period  T.  As  a  result  the  frequency  spectrum  of  the  OFDM  symbol  consists  of  the 
convolution  of  dirac  delta  functions  at  frequencies  corresponding  to  the  individual 
subcarrier  frequencies  with  the  frequency  spectrum  of  a  unit  amplitude  square  wave  with 
duration  T.  The  frequency  spectrum  of  a  unit  amplitude  square  wave  is  sinc(  njT),  which 
is  zero  at  all  frequencies  that  are  a  multiple  of  1/72  As  a  result  the  spectrum  of  an  OFDM 
symbol  with  overlapping  sine  functions  is  as  shown  in  Figure  2.4.  Note  that  at  the 
maximum  of  each  sine,  i.e.  frequency  spectrum  of  each  subcarrier,  all  other  sine 
functions  are  zero.  This  is  the  point  where  the  OFDM  receiver  calculates  the  spectral 
values  for  the  individual  subcarriers.  This  means  that  the  OFDM  receiver  is  able  to 
demodulate  each  subcarrier  free  of  interference  from  other  carriers.  Inter  Carrier 
Interference  (ICI)  would  exist  and  degrade  the  received  signal  if  this  were  not  the  case. 
Note  that  this  analysis  is  done  in  the  frequency  domain  of  the  OFDM  signal  and  therefore 
it  is  ICI  rather  than  Inter  Symbol  Interference  (ISI)  that  is  avoided.  [Ref.  4] 


Figure  2.4.  Subcarrier  Orthogonality  via  Subcarrier  Spectra  [after  Ref.  4]. 
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The  real  bandpass  signal  of  Eq.  2.2  can  also  be  expressed  in  complex  baseband 


form  as 


(2.3) 


The  complex  baseband  signal  is  the  inverse  Fourier  transform  of  N  QAM 
symbols.  Specifically  each  OFDM  symbol  is  the  inverse  Fourier  transform  of  the 
columns  or  blocks  of  Fig.  2.2.  In  practice  the  Inverse  Fast  Fourier  Transfonn  (IFFT)  is 
applied  to  the  QAM  symbols  to  generate  the  discrete  time  OFDM  signal  written  in 
equation  form  as 


(2.4) 


where  N  is  the  size  of  the  IFFT  and  is  equal  to  2 *  N  .[Ref.  2] 

Since  the  output  of  the  IFFT  must  produce  a  real  signal,  the  OFDM  blocks  of  Fig. 
2.2  must  be  modified  prior  to  application  of  the  IFFT.  From  Fourier  transfonn  theory  it  is 
known  that  a  complex  frequency  domain  signal  that  is  conjugate  symmetric  will  produce 
a  real  time  domain  signal  upon  application  of  the  inverse  Fourier  transform.  Therefore 
since  the  N  QAM  symbols  in  each  OFDM  block  represent  the  information  to  be 
transmitted,  and  are  defined  in  the  frequency  domain,  the  blocks  are  conjugated  and 
repeated  to  produce  the  conjugate  symmetric  counterpart  of  the  original  OFDM  block. 
The  result  is  an  OFDM  block  of  length  N(= 2*  N  )  complex  QAM  symbols.  The  first  N 
QAM  symbols  represent  the  real  frequency  portion  of  the  frequency  spectrum  and  the 
remaining  N  QAM  symbols  represent  the  imaginary  frequencies  in  the  spectrum.  Then 
upon  application  of  the  IFFT  to  each  block  of  length  N  complex  QAM  symbols,  a  block 
of  length  N  real  time  domain  samples  is  generated  which  are  converted  to  serial  in 
preparation  for  transmission. 
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2. 


Guard  Time  and  Cyclic  Extension 


One  of  the  major  advantages  of  OFDM  is  its  robustness  to  multipath  propagation 
and  the  associated  delay  spread  that  occurs.  The  delay  spread  is  the  time  difference  of 
arrival  of  transmitted  signals  that  follow  different  paths  from  the  source  to  the  receiver 
and  therefore  have  different  propagation  times. 

Dividing  the  input  data  stream  into  N  subcarriers  reduces  the  symbol  duration  by 

1/ N  .To  combat  the  effects  of  ISI  a  guard  time  is  incorporated  into  each  OFDM  symbol. 
The  duration  of  the  guard  time  is  chosen  to  be  at  least  as  long  as  the  maximum  expected 
delay  spread.  The  maximum  expected  delay  spread  is  the  time  difference  of  arrival 
between  the  first  and  last  transmitted  signals  to  arrive  at  the  receiver.  This  ensures  that 
multipath  receptions  of  one  symbol  do  not  interfere  with  the  receptions  of  the  following 
symbol. 

One  possibility  for  the  guard  time  is  to  transmit  no  signal.  However  this  would 
result  in  ICI,  which  is  crosstalk  between  carriers  that  results  from  a  loss  of  orthogonality. 
Recall  that  orthogonality  means  that  all  subcarriers  have  an  integer  number  of  cycles  in 
the  FFT  period  with  adjacent  subcarriers  having  a  number  of  cycles  that  differ  by  one. 
Integration  of  a  subcarrier  which  has  a  non-integer  number  of  cycles  in  the  FFT  interval 
will  result  in  a  loss  of  orthogonality  and  therefore  signal  degradation  due  to  ICI.  This 
effect  is  illustrated  in  Figure  2.5  below. 

Figure  2.5  shows  2  subcarriers  with  the  second  subcarrier  delayed.  Note  that  the 
guard  time  consists  of  a  zero  signal.  Also  note  that  subcarrier  one  has  exactly  two  cycles 
within  the  FFT  interval  while  subcarrier  two,  due  to  its  delay,  does  not  contain  an  integer 
number  of  cycles  within  the  FFT  interval  which  will  cause  ICI  as  discussed  above  and 
indicated  in  Fig.  2.5.  The  shaded  portion  of  subcarrier  2  indicates  the  region  of  the  signal 
that  causes  ICI.  One  can  see  that  if  the  subcarriers  were  extended  cyclically  in  the  guard 
interval  that  there  would  be  an  integer  number  of  cycles  in  the  FFT  interval  regardless  of 
the  multipath  delay  of  subcarrier  2.  The  cyclic  extension  is  implemented  by  simply 
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extending  the  length  of  each  OFDM  block  following  the  IFFT  operation  by  an  amount 
greater  than  or  equal  to  the  maximum  expected  delay  spread. 


Figure  2.6  illustrates  the  cyclic  extension  of  the  subcarriers  in  the  guard  interval. 
The  benefit  of  using  a  cyclic  extension  in  the  guard  interval  is  shown  in  Figure  2.7,  which 
shows  a  2  ray  multipath  environment  with  the  amplitude  of  the  second  ray  half  the 
amplitude  of  the  first  ray.  The  second  ray  arrives  after  the  first  ray  but  its  delay  spread  is 
less  than  the  duration  of  the  guard  interval  thus  there  is  no  1ST 
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[after  Ref.  4]. 


In  addition,  it  can  be  seen  that  there  are  an  integer  number  of  cycles  for  rays 
within  the  FFT  integration  time  which  maintains  the  orthogonality  of  the  subcarriers  and 
prevents  ICI.  The  subcarriers  in  Fig.  2.7  are  BPSK  modulated.  Therefore  only  two 
possible  phases  are  possible  and  they  differ  by  180  degrees.  The  phase  shifts  that  result 
are  indicated  in  the  figure.  Note  that  the  OFDM  receiver  would  actually  see  the 
summation  of  all  the  received  paths  but  the  subcarriers  are  shown  separately  for 
illustration  purposes. 
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environment  [after  Ref.  4], 


As  discussed  above,  the  maximum  delay  spread  must  not  exceed  the  guard 
interval  incorporated  in  the  OFDM  signal  or  ISI  and  a  loss  of  orthogonality  will  result 
causing  ICI.  Van  Nee  and  Prasad  [Ref.  4]  show  that  with  the  delay  spread  greater  than  the 
guard  interval  by  less  than  about  three  percent,  the  signal  interference  is  acceptable  but  at 
around  10%  delay  spread  excess  the  interference  is  unacceptable  for  communication 
purposes.  However  OFDM  is  still  significantly  more  robust  to  delay  spreads  that  exceed 
the  value  for  which  the  system  was  designed  relative  to  that  for  single  carrier  systems  as 
error  propagation  causes  abrupt  performance  degradation  in  single  carrier  systems  [Ref. 

4]- 


14 


B. 


CHANNEL  ANALYSIS 


In  general  communication  systems,  the  input  is  considered  to  be  the  transmitted 
signal  and  the  output  is  the  received  signal.  The  received  signal  is  the  transmitted  signal 
modified  by  the  channel  and  corrupted  by  additive  noise.  The  channel  is  usually  band 
limited  and  assumed  to  be  linear.  It  may  be  either  time  varying  or  time  invariant.  This 
concept  is  illustrated  in  Figure  2.8  below,  where  x(n)  is  the  input  signal,  y(n)  is  the 
output  signal,  h(n)  is  the  channel  impulse  response  and  w(n)  is  additive  gaussian  noise. 


In  OFDM,  each  of  the  N  subchannels  is  transmitted  with  its  own  carrier.  Consider 
an  OFDM  symbol  defined  as 


(2.5) 


where  k  refers  to  the  OFDM  block  or  symbol.  As  discussed  above,  must  be  modified 
to  be  conjugate  symmetric  prior  to  applying  the  IFFT.  Therefore  becomes 
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(2.6) 


where  N  =  2*  N  .  The  IFFT  then  modulates  the  N  subcarriers  resulting  in  the  real  time 
domain  signal  Xk , 


(2.7) 


The  channel  is  assumed  to  have  a  finite  impulse  response  of  length  L.  Therefore  a 
cyclic  extension  of  length  L  is  applied  to  the  OFDM  symbol  x*.  The  OFDM  signal 
becomes 


(*)  ik\  x,(k)---xN(k)] 


~^r 


cyclic  extension  IFFT^j,} 


Transmitted  OFDM  Signal 


(2.8) 


The  output  from  the  channel  is  the  convolution  of  the  channel  impulse  response 
with  the  input  signal.  Written  as  an  equation  this  is 
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(2.9) 


T/J 

In) 

=  A7<  < 

In)®  h(n) 

where  ®  represents  circular  convolution  and  yk(n)  is  the  kth  received  OFDM  block.  It  is 
well  known  that  for  continuous  time  signals  the  linear  convolution  in  the  time  domain  is 
equivalent  to  multiplication  in  the  frequency  domain.  In  order  for  this  to  hold  in  discrete 
time  one  of  two  conditions  must  be  met.  Either  the  block  length  N  must  be  infinite  or  at 
least  one  of  the  convolved  sequences  must  be  periodic  with  period  N.  Observe  that  the 
input  OFDM  signal  is  periodic  with  period  N  due  to  the  application  of  the  cyclic 
extension.  Therefore  the  received  frequency  domain  signal  is  the  product  of  the  input 
OFDM  frequency  domain  signal  and  the  frequency  response  of  the  channel.  In  equation 
form  the  kth  received  OFDM  symbol,  with  the  subscript  k  omitted  for  simplicity,  is 


Y.=HX 


i  i 


for  i  =  1,2. 


(2.10) 


where  the  subscript  i  refers  to  the  subcarrier.  Note  that  the  use  of  the  cyclic  prefix  reduces 
the  achievable  data  rate  by  a  factor  of  +  ^  .  However  for  most  applications  N»L  so 

the  loss  is  negligible  relative  to  the  identified  benefits  of  using  the  cyclic  prefix.  Most 
OFDM  systems  use  the  cyclic  prefix  [Ref.  5  and  6]. 


The  received  signal  is 


h(*)....*J,(t)]0[*(l).-*(L)]+w 


Received  OFDM  Signal 


(2.11) 
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As  indicated  in  Eq.  2.1 1  the  last  N  values  of  Vk  are  the  convolution  of  the  OFDM 
input  signal  with  the  channel  impulse  response.  Therefore  the  received  QAM  symbols  in 
the  kth  block  are  written  as  Eq.  2.12  where  the  last  N  QAM  symbols  are  the  complex 
conjugate  of  the  first  N  QAM  symbols  and  are  discarded  in  the  receiver. 


(2.12) 

Therefore,  the  received  QAM  symbol  for  the  ith  subchannel  and  the  k,h  block  can 
be  expressed  as 


where  H  is  defined  as 


(2.14) 

and  h(n)  is  zero  padded  to  length  N  in  taking  the  IFFT.  As  a  result  the  effect  on  the 
channel  is  that  it  is  also  divided  into  narrow  subchannels  corresponding  to  the  subcarriers. 

For  large  N  ,  the  subchannels  are  very  narrow  and  there  is  negligible  change  in  H t  and 
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versus  frequency.  In  other  words  the  frequency  selective  fading  that  may  be 

present  is  eliminated.  Subchannel  independence  can  be  shown  provided  the  additive  noise 
is  Gaussian( whether  the  noise  is  colored  or  not)  [Ref.  3].  This  result  is  illustrated  in  Fig. 
2.9  below. 

The  rectangles  represent  the  subchannels  that  are  used  by  the  subcarrier  bands. 
Each  subchannel  is  assumed  flat  for  large  N.  In  Fig.  2.9,  only  1 1  subcarriers  exist 
however  in  practice  the  number  of  subcarriers  is  typically  between  500  and  2000.  With 
this  in  mind  one  can  see  that  the  frequency  response  within  each  of  the  subchannels 
would  be  flat  as  the  width  of  the  subchannel  would  be  drastically  reduced. 


19 


THIS  PAGE  INTENTIONALLY  LEFT  BLANK 


20 


III.  OFDM  SIMULATION  MODEL 


The  development  of  a  computer-based  simulation  of  an  underwater  OFDM 
communication  system  is  the  primary  thrust  of  this  thesis.  The  model  simulates  most  of 
the  significant  hardware  that  would  exist  in  an  actual  system.  This  chapter  describes  the 
details  of  the  model  and  the  theory  behind  the  methods  applied.  The  intent  is  not  to 
exhaustively  model  every  component  in  an  actual  system  but  to  model  those  components 
and  features  of  a  real  system  that  contribute  the  dominating  characteristics  of  the 
system’s  performance.  It  is  also  understood  to  be  unfeasible  and  foolish  to  attempt  to 
develop  a  model  from  scratch  that  will  have  the  ability  to  experiment  with  every  possible 
design  and  operating  feature  of  the  system.  With  that  in  mind,  the  model  was  developed 
to  be  comprehensive  enough  to  explore  some  of  the  basic  specifications  common  to  any 
communication  system  immediately  but  then  to  also  provide  a  test  bed  upon  which 
further  modeling  and  experimentation  could  take  place. 


A.  METHODOLOGY 

The  development  of  each  block  in  the  model  followed  essentially  the  same  path. 
First  the  topic  of  interest  was  explored  in  the  literature  to  the  extent  required.  The  extent 
of  literature  research  was  a  function  of  the  number  of  choices  available  to  accomplish  the 
task,  the  experience  of  this  author  in  the  subject  and  the  relative  importance  of  the  block 
to  the  entire  model.  For  example,  the  modulation  and  demodulation  blocks  were  only 
mildly  researched  while  the  topic  of  forward  error  correction  was  researched  much  more 
extensively.  Figure  3.1  provides  a  block  diagram  of  the  entire  simulation. 
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B. 


TRANSMITTER 


The  transmitter  covers  all  aspects  of  signal  processing  and  transmission  from 
information  signal  development  to  Quadrature  Amplitude  Modulation  (QAM),  to  Double 
Side  Band  Suppressed  Carrier  (DSBSC)  modulation  of  the  signal.  Each  of  the  transmitter 
blocks  of  Fig.  3.1  is  discussed  in  the  following  numbered  subsections. 


1.  OFDM  Parameter  Determination 


In  designing  a  communication  system  there  are  several  parameters  of  interest  that 
must  be  determined.  The  final  design  will  be  a  compromise  between  conflicting 
requirements  to  achieve  the  optimal  perfonnance.  Parameter  determination  is  performed 
in  the  simulation  by  the  Matlab  function  parameters. m  contained  in  Appendix  C.  In 
general  there  are  three  main  parameters  to  determine:  bit  rate,  R ,  channel  delay  spread, 
tds,  and  bandwidth,  W.  Delay  spread  is  a  function  of  the  channel  and  therefore  represents 
the  opening  argument  in  designing  the  system.  The  OFDM  system  must  be  designed  such 
that  the  tolerable  delay  spread  of  the  system  is  greater  than  the  delay  spread  of  the 
channel  to  avoid  signal  degradation.  Therefore  the  guard  time,  tg,  is  typically  two  to  four 
times  the  maximum  expected  delay  spread.  That  is, 


(3.1) 


The  delay  spread  for  the  underwater  acoustic  sound  channel  is  a  function  of  but 
not  limited  to  bottom  type,  transmission  range,  transmitter  and  receiver  depths  and 
frequency  band.  Coatelan  and  Glavieux  [Ref.  1]  report  measured  delay  spreads  in  Table 
3.1  below.  The  table  identifies  the  transmission  range  between  the  receiver  and  the 
transmitter,  the  seafloor  composition  and  the  associated  tds  measured. 
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Bottom 

Transmission 

Delay  Spread 

Type 

Range  (m) 

(msec) 

Rocky 

<  1000 

<  10 

Sludgy 

1500-2500 

<10 

Sandy 

1500-2500 

>20 

Unknown 

1430 

~20 

Table  3.1.  Typical  Delay  Spread  Values. 


The  next  step  is  to  determine  the  OFDM  symbol  duration.  The  OFDM  symbol 
duration  is  the  sum  of  the  cyclic  prefix  duration  and  the  FFT  interval  duration.  Since  the 
signal  in  the  cyclic  prefix  does  not  contain  any  new  information  it  is  desirable  to  have  the 
FFT  interval  duration  much  greater  than  the  guard  time  to  maximize  the  achievable  bit 

rate.  However  there  is  a  limit.  As  T  increases,  the  subcarrier  spacing,  A fsc  ?  decreases  for 

a  given  bandwidth.  Smaller  subcarrier  spacing  results  in  greater  sensitivity  to  phase  errors 
and  frequency  offset.  In  practice  the  symbol  duration  is  typically  5  or  more  times  tg, 
therefore 


T  >  : 

5  *  t 

g 

For  channels  with  frequency  selective  fading,  smaller  A fsc  values  are  preferred 
as  this  ensures  the  subchannel  frequency  response  is  flat.  Conversely,  for  channels  with 
large  Doppler  spreading  characteristics  larger  A fsc  values  are  better. 

Figure  3.2  below  presents  correlator  output  results  measured  by  Coatelan  and 
Glavieux  [Ref.  1].  The  results  in  Fig.  3.2  correspond  to  the  1430  m  acoustic  channel  of 
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Table  3.1  with  a  sludgy  seafloor  and  a  depth  of  90  m.  Note  that  the  delay  spread  in  Fig. 
3.2  is  approximately  20  msec  and  the  first  path  received  is  not  necessarily  the  most 
energetic  one  [Ref.  1]. 


With  the  guard  time  and  symbol  duration  established,  the  number  of  subcarriers 
N  can  be  determined.  There  are  two  approaches  that  can  be  taken.  The  first  is  starting 
with  a  known  bandwidth  constraint.  In  this  case  the  number  of  subcarriers  is  simply  the 
available  bandwidth,  W,  divided  by  the  subcarrier  spacing,  therefore 


A 4  [  (33) 

The  other  approach,  which  is  the  one  used  in  the  simulation,  is  to  start  with  a 
desired  bit  rate,  R.  Then  the  number  of  subcarriers  is  simply  R  divided  by  the  number  of 
bits  transmitted  by  each  subcarrier,  bsc.  That  is, 
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(3.4) 


a.  Example  Calculation 


Suppose  that  a  system  is  to  be  designed  that  would  provide  a  10  kbps  data 
rate  and  the  maximum  expected  delay  spread  for  the  acoustic  channel  of  interest  is  27 
msec.  That  is, 

R  =  10  kbps 
tds  =  27  msec 

let  tg  =  4 tds  =  4(27)  =  108  msec 
let  T  =  6tg  =  6(l 08)  =  648  msec 


Assume  that  the  system  is  capable  of  transmitting  6  bits  per  subcarrier, 
which  corresponds  to  64-QAM,  then 


N  = 


10  £3 
6 


=  1667 


subcarriers 


However  since  the  QAM  symbols  on  each  of  the  subcarriers  will  be 
modulated  using  the  IFFT,  it  is  computationally  more  efficient  to  have  the  number  of 

subcarriers  equal  to  a  power  of  two.  Therefore  N  is  set  to  2048  subcarriers  but  only  1667 
of  the  subcarriers  have  non-zero  values  with  the  remaining  subcarriers  providing 
oversampling  to  avoid  aliasing.  In  essence  the  signal  itself  is  unaffected  and  changing  N 
only  affects  the  internal  signal  processing  efficiency.  Next  the  FFT  interval  duration  is 
calculated  as, 


TpfT  —  T  —  tg  —  648  — 108  —  540  msec 
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Another  requirement  is  that  there  be  an  integer  number  of  samples  within 
both  the  FFT  interval,  Tfft,  and  the  entire  OFDM  symbol,  T,  in  order  to  maintain 
orthogonality  of  the  subcarriers.  For  this  example  there  must  be  2048  samples  within  the 
FFT  interval.  The  FFT  sampling  rate,  Fs_Fft,  is  then, 


F 


N  2048 


s.FFT 


1  FFT 


540 


=  3.793 


However  the  number  of  samples  within  the  OFDM  symbol  interval,  T,  is, 

nT  =  Fs  FFT  *  T  =  3.793  *  540  =  2457.6 


which  is  not  an  integer.  Therefore  the  subcarriers  will  not  be  orthogonal  and  signal 
degradation  will  occur.  The  remedy  is  to  adjust  the  sampling  rate  slightly.  The  sampling 
rate  is  increased  to  achieve  an  integer  number  of  samples  within  the  FFT  interval  and  the 
OFDM  symbol  interval.  For  this  example  FSyFFT becomes  3842.6  Hz,  which  yields 

nT  =  Fs  FFT  *  T  =  3842.6  *  648it  —  3  =  2490  samples 
nFFT  =  Fs  fft  *  Tfft  =  3842.6  *  540if  —  3  =  2075  samples 
ng  =  FsFFT  =  3842.6 *  108ii  —3  =  415  samples 


where  nT,nFFT  and  ng  are  the  number  of  samples  in  the  OFDM  symbol  interval,  the  FFT 

interval  and  the  guard  interval  respectively.  Note  that  the  OFDM  symbol  is  oversampled 
by  increasing  the  sampling  rate  since  there  are  2075  vice  2048  samples  within  the  FFT 
interval.  Synchronization,  which  will  be  discussed  in  the  receiver  section,  will  detennine 
which  of  the  2075  samples  in  the  FFT  interval  are  to  be  used  for  signal  reconstruction.  As 
a  result  of  changing  the  sampling  rate, 


T  =  =  0.5330  sec 

f1 

1  s,FFT 


tg  =  — —  =  0.1 15  sec 

Fs,FFT 
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and  finally, 


AT*  =7^  =1.876  Hz 
W  =  N  *  Afsc  =2.0264  kHz 

2.  Signal  Generation 


In  underwater  communication  several  possible  data  forms  may  be  transmitted 
including  video,  voice  and  text  fonnat  messages.  For  the  simulation,  voice  data  in  .wav 
format  is  used  as  the  information  signal.  Voice  data  was  chosen  as  it  requires  an  analog 
signal  to  be  converted  to  binary  format  via  sampling,  quantization  and  decimal  to  binary 
conversion.  Once  the  data  is  in  binary  form  the  original  form  of  the  data  is  irrelevant  to 
the  OFDM  system.  Figure  3.3  presents  the  information  signal  used  in  the  simulation.  The 
upper  plot  is  the  discrete  time  signal.  The  sampling  rate  of  the  analog  voice  signal  is 
1 1 .025  kHz  and  the  data  is  quantized  using  8  bits  per  sample.  The  lower  plot  is  the 
frequency  spectrum  of  the  transmitted  signal. 


Bin  number  (k) 


Figure  3.3.  Information  Signal. 
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3. 


Quantization 


The  sampled  information  signal  is  then  quantized  using  uniform  quantization.  The 
level  of  quantization  is  adjustable  but  8-bit  quantization  is  primarily  used  in  the 
simulation  which  results  in  256  quantization  levels.  Quantization  is  perfonned  by  the 
Matlab  function  quantization.m  contained  in  Appendix  C. 

The  output  of  the  quantizer  is  a  stream  of  decimals  from  0  to  M- 1  where  M  is  2blts 
and  bits  is  the  number  of  bits  per  binary  digits  per  quantization  level.  The  decimal  value 
corresponds  to  the  quantization  level  for  which  the  sampled  value  of  the  information 
signal  was  nearest.  The  quantizer  also  outputs  a  codebook  that  identifies  the  quantization 
level  corresponding  to  the  decimal  value  for  signal  reconstruction.  Figure  3.4  clarifies  the 
quantization  process  with  a  simple  2-bit  example. 


2  Bit 

4  Level 

+  00 

A 

i 

-  3  — 

-  X.XM 

Partition  3 

9 

x.xxx 

Partition  2 

vwvwv 

-  1  — 

- 

Partition  1 

-  0  — 

-  X.ffi 

1 

r 

Decimal  Quantization 

00 

Value 

Level 

F igure  3.4.  Quantization  Illustration. 
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4. 


Decimal  to  Binary  Conversion 


Decimal  to  binary  conversion  is  performed  by  the  Matlab  code 
binaryconversion.m  contained  in  Appendix  C.  The  process  involves  straightforward 
conversion  of  a  base  10  integer  to  the  corresponding  base  2  representation. 

5.  QAM  Buffer 

QAM  modulation  involves  the  mapping  of  bits  to  a  QAM  constellation.  The  size 
of  the  constellation,  Q,  is  a  function  of  the  number  of  bits  per  QAM  symbol. 

Q  =  2  7  (3.5) 

In  general  the  number  of  bits  per  QAM  symbol  will  not  be  the  same  as  the 
number  of  bits  per  quantization  level.  Therefore  it  is  necessary  to  buffer  the  bits  prior  to 
mapping  to  a  QAM  constellation.  The  Matlab  code  qam  buffer.m  contained  in  Appendix 
C  performs  the  buffering  operation.  This  involves  converting  the  bit  stream  from  the 
decimal  to  binary  conversion  process  into  a  matrix  of  bits.  Each  row  in  the  matrix 
represents  a  QAM  symbol  in  binary  fonn  and  the  number  of  columns  per  row  is  q.  The 
number  of  rows  is  a  function  of  the  size  of  the  data  to  be  transmitted.  Since  the  matrix 
must  be  square  in  order  to  convert  all  bits  to  a  QAM  symbol,  up  to  q- 1  zeroes  are  added 
to  the  data. 


6.  Reed  Solomon  Encoder 

Due  to  the  multipath,  frequency  selective  fading  nature  of  the  underwater  acoustic 
sound  channel  the  subcarriers  will  arrive  with  different  amplitudes.  Some  subcarriers  may 
experience  severe  enough  fading  that  they  are  lost  entirely.  As  a  result  the  bit  error  for 
those  subchannels  will  approach  !4  and  will  dominate  the  BER  performance  of  the  entire 
system[Ref.  4], 

Forward  error  correction  coding  is  used  to  combat  the  severe  fading  affect.  The 
coding  is  applied  across  all  subcarriers  such  that  the  errors  of  the  weak  subcarriers  can  be 
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corrected.  The  performance  of  the  system  is  then  a  function  of  the  average  received 
power  vice  the  received  power  of  the  weak  subcarriers  [Ref.  2]. 

In  the  simulation  block  coding  is  applied,  specifically  Reed  Solomon  block 
coding  which  is  perfonned  at  the  symbol  level.  Coding  involves  the  addition  of  redundant 
symbols  to  increase  the  Hamming  distance  between  symbols.  In  general  a  message  of 
length  k  is  added  thereby  creating  a  code  word  of  length  n.  The  coding  rate  r  is  defined  as 
k/n.  The  minimum  Hamming  distance,  dmuh  is  the  minimum  number  of  different  symbols 
between  any  2  code  words. 

For  a  given  code  with  an  associated  dmin,  t  errors  can  be  corrected  by  the  code 
where  t  is 


where  the  floor  function  denotes  rounding  down  to  the  nearest  integer[Ref.  2].  The 
minimum  Hamming  distance  has  an  upper  bound  defined  by  the  number  of  redundant 
symbols  as 

(3.7) 

In  terms  of  bit  level  codes,  only  repetition  and  single  parity  check  codes  reach  this 
upper  bound.  Reed  Solomon  is  a  class  of  non-binary  codes  that  reach  this  upper 
bound. [Ref.  4] 

Reed  Solomon  codes  are  defined  for  code  words  of  length  n  symbols  associated 
with  m  bits  per  symbol  where 

(3.8) 

The  required  message  length,  k,  is  related  to  message  and  the  number  of  bits  per 
symbol  as 

(3.9) 
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Associated  with  the  number  of  symbol  errors  the  code  can  correct  is  the  number 
of  bits  the  code  is  capable  of  correcting,  4,  defined  as 


<  m  *  floor 


r  n  - 


(3.10) 


Equation  3.10  is  true  provided  the  bit  errors  occur  within  the  maximum 
correctable  number  of  symbol  errors.  For  example  consider  a  code  designed  to  be  able  to 
correct  3  symbol  errors  with  4  bits  per  symbol.  The  code  will  not  be  able  to  correct  an 
arbitrary  occurrence  of  4  bit  errors  as  they  may  occur  in  4  different  symbols.  As  a  result 
Reed  Solomon  codes  are  useful  for  ‘bursty’  error  channels  such  as  the  multipath, 
frequency  selective  fading  underwater  acoustic  channel  since  the  errors  are  concentrated 
in  a  few  of  the  subcarriers  that  are  deeply  faded. 


The  Reed  Solomon  encoder  receives  QAM  coded  decimals  in  serial  format  from 
the  QAM  buffer.  The  vector  of  decimal  symbols  must  be  reordered  to  a  k  column  square 
matrix.  Therefore  subsymbols  may  need  to  be  added  to  create  a  square  matrix,  in  which 
case  they  are  added  to  the  beginning  of  the  vector. 


7.  Block  Interleaver 

The  frequency  selective  fading  causes  varying  amplitudes  of  the  subcarriers 
where  some  of  the  carriers  will  be  unreliable  due  to  severe  fading.  As  a  result  burst  errors 
vice  random  errors  will  occur  over  the  weak  subcarriers.  Most  FEC  codes  are  designed 
for  random  errors  not  burst  errors[Ref.  4].  The  transmitter  pennutes  the  data  such  that  the 
adjacent  symbols  are  separated  by  several  blocks  after  interleaving.  The  receiver 
performs  the  inverse  operation. 

As  discussed  Reed  Solomon  codes  are  designed  to  correct  a  certain  number  of 
errors  per  block.  The  interleaving  must  distribute  the  subsymbols  across  several  blocks 
such  that  the  bursts  caused  by  the  severely  faded  subcarriers  appear  random  in  the 
receiver  after  de-interleaving.  Therefore  rather  than  having  bursts  of  symbol  errors  within 
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each  block  corresponding  to  the  severely  faded  subcarriers,  the  errors  within  each  block 
will  be  randomly  distributed  which  is  optimal  for  the  Reed  Solomon  decoding  algorithm. 

The  interleaving  process  is  essentially  a  reordering  of  the  data  and  is  performed 
by  the  reshape. m  command  in  Matlab. 

The  input  data  to  the  interleaver  is  an  n  column  matrix  where  the  data  is  ordered 
in  time  sequence  across  the  rows.  Therefore  the  sequential  order  of  the  data  in  the  matrix 
begins  at  row  1,  column  1,  then  proceeds  across  row  1  to  column  n,  then  to  row  2, 
column  1  etc. 

The  reshape  function  takes  the  data  column  wise  such  that  the  subsymbols  are 
interleaved  in  frequency.  The  process  is  illustrated  in  Fig.  3.5  with  n  equal  to  7  and  k 
equal  to  3.  In  the  simulation  the  length  of  the  input  vector  A  is  on  the  order  104,  n  is  at 
most  64,  and  the  number  of  active  subcarriers  per  OFDM  symbol  is  typically  less  than 
10  .  Therefore  the  interleaving  is  over  several  blocks  as  desired. 
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Figure  3.5.  Reed  Solomon  Encoder  and  Block  Interleaver. 


8.  QAM  Modulator 

QAM  Modulation  maps  the  QAM  coded  decimals  representing  points  on  a  QAM 
constellation  into  real  and  imaginary  components.  The  real  component  is  called  the  in- 
phase  component,  and  the  imaginary  component  is  called  the  Quadrature  component  of 
the  signal.  The  mapping  from  decimal  to  real  and  imaginary  components  is  illustrated  in 
Fig.  3.6  below. 
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Figure  3.6.  16  QAM  Constellation  with  In-phase  and  Quadrature  Mapping. 


The  Cartesian  coordinates  of  in-phase  and  Quadrature  components  are  converted 
to  polar  coordinates  prior  to  creating  the  OFDM  symbol  blocks.  The  QAM  Modulator 
creates  the  frequency  spectrum  of  the  signal  using  the  complex  QAM  symbols  now  in 
magnitude  and  phase  components.  The  goal  is  to  generate  blocks  of  complex  QAM 
symbols  such  that  applying  them  to  the  IFFT  operation  results  in  a  real  signal.  The  IFFT 
is  given  by 


(3.11) 


where  N  is  the  length  of  the  DFT.  The  discrete  time  signal  x„  will  be  real  provided  the 
following  condition  of  conjugate  symmetry  holds  for  the  frequency  domain  signal  A„ 
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(3.12) 


The  OFDM  blocks  are  then  fonned  by  sequentially  filling  the  active  subcarriers  in 
each  block  with  the  complex  QAM  symbols. 

In  the  parameter  determination  algorithm  the  number  of  active  subcarriers,  Naci , 

and  the  number  of  zero  subcarriers,  N:er ,  are  determined.  The  zero  subcarriers  are 
required  to  prevent  aliasing  from  signal  power  near  the  Nyquist  sampling  rate  which 

corresponds  to  N  .  The  number  of  data  blocks,  AWa  required  to  transmit  a  given  amount 
of  data  is  given  by  the  number  of  complex  QAM  symbols  to  be  transmitted  divided  by 
the  number  of  active  subcarriers  rounded  up  to  the  nearest  integer.  As  a  result  of 
rounding  up  to  get  the  N block,  there  is  a  surplus  of  active  subcarriers  relative  to  the  number 
of  QAM  symbols  to  be  transmitted.  In  future  work  this  surplus  could  be  filled  with 
protocol  information.  The  surplus  subcarriers  are  filled  with  zero  symbols  in  the 
simulation. 


The  original  design  placed  the  required  number  of  zeroes  at  the  beginning  of  the 
first  block.  However,  this  was  found  to  be  a  poor  choice  as  it  results  in  an  uneven  power 
distribution  across  the  frequency  spectrum.  The  current  design  spreads  the  zero  symbols 
across  the  frequency  band  used  for  transmission. 


to 


The  number  and  location  of  zero  symbols  to  be  inserted  are  determined  according 


(3.13) 


where  N0  is  the  total  number  of  zero  symbols  to  be  inserted  and  N0AMi  is  the  total 

number  of  complex  QAM  symbols  to  be  transmitted.  In  general  N0  will  not  be  a  multiple 

of  Nbiock,  therefore  the  number  of  zero  symbols  inserted  will  not  be  the  same  for  each 
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OFDM  block.  The  following  method  is  used  to  determine  the  number  of  zeroes  inserted 
into  each  OFDM  block. 


r 


Z  col  =  fl°°r 


N 


0 


V  ^  block  J 


Z  col  ,1  ^0  (N  block  1 )*Zcol 


(3.14) 


where  Zcol ,  is  the  number  of  zero  symbols  inserted  into  the  first  OFDM  block  and  Zco/  is 

the  number  of  zero  symbols  inserted  into  the  remainder  of  the  OFDM  blocks.  The 
process  is  illustrated  in  Fig.  3.7  below. 

Application  of  Eq.  3.12  to  the  N  QAM  symbols  in  each  of  the  OFDM  blocks 
results  in  OFDM  symbols  which  are  conjugate  symmetric.  The  conjugate  symmetry  will 
ensure  that  each  of  the  OFDM  blocks  of  complex  QAM  symbols  will  result  in  N(  =  2 N  ) 
real  samples  upon  application  of  the  IFFT.  The  result  is  shown  in  Fig.  3.8  with  the  same 

parameters  as  introduced  in  Fig.  3.7.  The  symbol  ^  in  Fig.  3.8  refers  to  complex 
conjugation. 
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Figure  3.7.  QAM  Modulator. 


The  QAM  Modulator  also  creates  the  pilot  blocks  which  are  used  for  channel 
estimation  and  equalization.  Although  the  theory  behind  channel  estimation  and 
equalization  will  be  discussed  later,  the  creation  of  the  pilot  blocks  is  presented  here. 

In  simple  terms  the  pilot  blocks  are  symbols  known  to  the  receiver  that  are  used  to 
estimate  the  channel  frequency  response.  In  the  simulation  two  pilot  blocks  are 
implemented  and  make  up  the  first  two  OFDM  symbols  that  are  transmitted.  The  pilot 
blocks  are  made  up  of  pseudo  random  complex  OFDM  symbols  which  are  detennined  as 
follows. 

The  data  used  to  create  the  pilot  blocks  comes  from  the  data  blocks  illustrated  in 
Fig.  3.7.  Each  complex  QAM  symbol  is  drawn  from  the  data  blocks  and  placed  in  the 
first  pilot  block  as  follows, 
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x P,Xc,r)=  Xh(r,c)  for 


r  —  N  act  1 5  N  act  ,  .  .  •  1 

< 

c  =  N  N  1  ( 

iN  block  J  iN  block  ?  *  *  *  1 

3.15) 

Equation  3.12  is  applied  to  the  first  N  complex  QAM  symbols  creating  a  length 
A  block  of  complex  QAM  symbols  that  is  conjugate  symmetric.  Reversing  the  order  of 
the  first  pilot  block  creates  the  second  pilot  block. 

The  final  effect  of  the  QAM  Modulator  is  that  by  adding  the  complex  conjugate 
counterpart  to  each  of  the  N  complex  QAM  symbol  blocks  the  input  data  is  effectively 
interpolated.  Typically  interpolation  involves  the  insertion  of  one  or  more  zeroes  in  a 
time  sequence  with  the  effect  of  increasing  the  sampling  rate  and  ‘compressing’  the 
frequency  spectrum  of  the  sequence.  In  the  QAM  Modulator  the  interpolation  constant  is 
2  as  the  number  of  data  points  is  doubled.  The  result  is  that  the  sampling  rate,  FSiFFT,  is 
doubled.  Therefore  the  new  sampling  rate  in  the  FFT  interval  is, 

(3.16) 
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9. 


IFFT/OFDM  Modulator 


The  IFFT  is  used  to  modulate  the  complex  QAM  symbols  to  the  subcarrier 
frequencies.  In  other  words  each  subcarrier  is  modulated  by  the  corresponding  complex 
QAM  symbol. 


The  IFFT  can  also  be  expressed  in  matrix  form  as 


(3.17) 


(3.18) 


(3.19) 


where  W_N  e  CNxN  and  X_N  e  C  VxM  therefore  xN  e  RNxM  for  X_N  conjugate  symmetric 
as  defined  in  Eq.  3.12. 


Note  that  C  refers  to  the  complex  subspace  and  R  refers  to  the  real  subspace  and 
the  M  refers  to  the  number  of  columns  in  the  matrix.  In  the  simulation  M  is  Nbiock ■  The 
output  from  the  IFFT  is  a  matrix  of  real  samples  with  each  column  representing  a 
baseband  real  OFDM  symbol.  Columns  1  and  2  are  the  pilot  symbols  and  the  remaining 
columns  contain  the  infonnation  to  be  transmitted. 
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10.  Cyclic  Extension 


The  need  for  the  cyclic  extension  was  discussed  in  Chapter  II.  The 
implementation  of  the  cyclic  prefix  in  the  simulation  is  quite  simple.  First  the  required 
length  of  the  cyclic  extension  is  determined  by 

L  =  round (tg  (3  20) 

where  L  is  the  number  of  discrete  time  samples  that  must  be  prefixed  to  each  real  OFDM 
symbol. 

Recall  that  xN  from  the  IFFT  is  an  N  x  N block  matrix  of  real  time  samples 

therefore  the  cyclic  extension  results  in  a  matrix  that  is  N+L  x  Nbiock  with  the  first  L 
samples  being  identical  to  the  last  L  samples  in  each  column.  Therefore  the  resulting 
matrix  is 


The  columns  of  the  matrix  of  Eq.  3.21  represent  the  real  baseband  OFDM  signals 
that  are  to  be  transmitted. 
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11.  Parallel  to  Serial  Converter 


The  parallel  to  serial  converter  simply  reshapes  xNc  to  a  vector  of  length 

(. Nbiock+2)*(N+L ).  In  the  simulation,  round-off  error  causes  the  output  from  the  IFFT  to 
contain  small  imaginary  components.  Therefore  the  real  part  of  the  serial  time  samples  is 
taken. 


12.  Double  Side  Band  Modulator 


After  serial  conversion  the  signal  is  baseband  and  must  be  modulated  to  the 
desired  transmission  band.  The  modulation  is  performed  in  discrete  time  to  take 
advantage  of  the  ability  to  use  high  order  digital  filters. 


In  order  to  prevent  aliasing  the  signal  must  be  interpolated.  The  interpolation 
factor,  /,  is  determined  as 


(3.22) 


where  fc  is  the  carrier  frequency. 


Interpolation  is  performed  by  inserting  I- 1  zeroes  between  each  sample.  The 
interpolated  signal  can  be  expressed  as 


where  x(n)  is  the  baseband  signal  and  y(m)  is  the  baseband  signal  following  interpolation 
by  the  factor  I.  The  increased  sampling  rate  is  then 


F 


s ,  mod 


(3.24) 
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Figure  3.9  illustrates  the  effect  of  interpolation  on  the  signal.  The  signal  of  Fig.  3.9  is  for 
a  5  kbps  transmission  rate  using  16-QAM.  There  are  1024  subcarriers  with  750  active 
subcarriers  which  results  in  a  bandwidth,  W,  of  1.501  kHz.  The  baseband  signal  is 
interpolated  by  a  factor  of  4  as  shown  in  the  lower  graph  of  Fig.  3.9. 


Baseband  OFDM  Signal 


Frequency  (Hz) 


Figure  3.9.  Interpolated  Baseband  OFDM  Signal. 


The  interpolated  signal,  y(m),  must  be  filtered  to  remove  the  unwanted  spectral 
energy  above  the  OFDM  signal  bandwidth,  W,  as  seen  in  Fig.  3.9.  A  digital  Finite 
Impulse  Response  (FIR)  filter  is  used.  The  FIR  filter  is  designed  using  a  Hamming 
window  with  the  following  parameters. 
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The  filter  parameters  are  shown  for  a  general  FIR  filter  designed  using  the 
Hamming  window  method  in  Fig.  3.10.  The  variable  co,  is  the  digital  frequency  defined 
as 


(3.25) 


where/ refers  to  the  analog  frequency  and  Fs  is  the  sampling  rate  used  to  sample  the 
continuous  time  signal.  The  output  of  the  FIR  filter  is  defined  as, 

yf{n)=h{n)*y{n)  ( 

where  h(n)  is  the  impulse  response  of  the  FIR  filter  and  symbol  *  refers  to  linear 
convolution. 


Figure  3.11  shows  the  effect  of  the  FIR  filter  on  the  interpolated  baseband 
spectrum  shown  in  the  lower  graph  of  Fig.  3.9.  The  interpolated  baseband  frequency 
spectrum  is  repeated  in  the  upper  graph  of  Fig.  3.11.  The  lower  graph  shows  the 
frequency  spectrum  of  the  interpolated  baseband  frequency  spectrum  following  filtering. 
The  transmission  parameters  are  the  same  as  for  Fig.  3.9.  The  frequency  response  of  the 
FIR  filter  is  shown  in  Fig.  3.12. 

After  interpolation  and  filtering,  the  signal  is  ready  for  modulation  to  the  desired 
carrier  frequency.  The  modulation  is  performed  by  discrete  time  domain  multiplication 
with  a  cosine  modulating  signal.  The  modulating  signal  is  defined  as 
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(3.27) 


COS 


2  7tfc  n  ^ 
F 

\  s  ,  mod  y 


where  n  is  defined  from  1  to  the  number  of  samples  contained  in  the  signal  to  be 
modulated,  y/n). 


Figure  3.10.  General  FIR  Filter  (Hamming  Window  Design). 


46 


The  frequency  spectrum  of  the  cosine  modulating  signal  consists  of  dirac  delta 
functions  at  the  absolute  value  of  the  carrier  frequency  on  both  the  real  and  the  imaginary 
frequency  axis.  The  modulation  is  expressed  in  discrete  time  as 


Tmod(«)=  m{n)xyf(n) 


(3.28) 


Figure  3.11.  Interpolated  Baseband  OFDM  Signal  following  FIR  Filter. 


The  effect  of  modulation  onyyis  shown  in  Fig.  3.13.  The  upper  graph  in  Fig.  3.13 
repeats  the  frequency  spectrum  of  the  baseband  OFDM  signal  prior  to  interpolation.  The 
center  graph  is  the  frequency  spectrum  of  the  modulating  signal  which  illustrates  the 
dirac  delta  functions  at  the  carrier  frequency.  Finally  the  lower  graph  is  the  frequency 
spectrum  of  the  modulated  OFDM  signal.  This  is  the  form  of  the  signal  that  is  transmitted 
through  the  channel.  In  practice  the  signal  would  be  applied  to  a  Digital  to  Analog 
converter  (DAC)  and  an  associated  low  pass  filter  (LPF).  For  the  simulation  these 
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components  are  assumed  to  be  a  part  of  the  channel  transfer  function.  Similarly,  an 
Analog  to  Digital  Converter  (ADC)  and  associated  LPF  would  be  used  in  practice  to 
return  the  signal  to  discrete  time  following  transmission  through  the  channel.  These  are 
also  assumed  to  be  included  in  the  channel  transfer  function. 
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Figure  3.12.  FIR  Filter  Frequency  Response. 


48 


m 

^3 


*5 

£ 

H 

Uh 

P-. 


03 

■T3 


100 

50 

0 

-50 


H 

Oh 


10 

-loh 


>5  -30  h 


-50 


Frequency  (Hz) 


Real  Frequencies 


2 


Imaginary  Frequencies 


A 


i 


0  2000  4000  6000  8000  10000  12000  14000  16000  18000 

Frequency  (Hz) 


2000  4000  6000  8000  10000  12000  14000  16000  18000 

Frequency  (Hz) 


Figure  3.13.  Modulation  Signals. 


The  DSB  Modulator  completes  the  transmitter  for  the  simulation.  The  next  effect 
on  the  signal  in  practice  would  be  the  channel.  However  in  the  simulation  the  channel 
transfer  function  is  not  applied  until  after  demodulation  in  the  receiver  for  reasons  that 
will  be  discussed  later. 

C.  MMPE  CHANNEL  MODEL 

MMPE  is  the  Monterey-Miami  Parabolic  Equation  model.  It  provides  an 
approximate  solution  to  the  Helmholtz  equation  based  on  the  efficient  split-step  Fourier 
(SSF)  marching  algorithm.  The  Helmholtz  equation  is  derived  from  the  solution  of  the 
linear,  three-dimensional,  lossless,  homogeneous  wave  equation  describing  the 
propagation  of  small  amplitude  acoustic  signals  in  an  ideal  fluid  [Ref.  6].  The  velocity 
potential  is  assumed  to  have  a  time-harmonic  solution,  or  in  other  words,  the  acoustic 
field  is  generated  by  a  continuous  wave  source.  The  value  of  a  time  harmonic  field  at  any 
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point  and  time  in  space  is  a  function  of  only  one  frequency.  Time  harmonic  solutions  are 
important  because  solutions  to  the  wave  equation  with  arbitrary  time  dependence  can  be 
found  through  application  of  Fourier  transform  techniques  to  the  time  hannonic  solutions 
[Ref.  7], 


The  following  development  follows  the  work  by  Houdeshell.  [Ref.  8]  The  linear 
wave  equation  is  given  in  Eq.  3.29  with  the  complex  acoustic  pressure  field  given  as  Eq. 
3.30, 


v 2  p  + 


l  d2p 

c 2  dt2 


0 


(3.29) 


(3.30) 


where  P  is  the  complex  acoustic  pressure  field  and  p  is  the  amplitude  of  the  acoustic 
pressure.  The  Helmholtz  equation  is  then 


(3.31) 


A  cylindrical  coordinate  system  is  used  in  the  solution.  Therefore  the  position 

vector,  Y  ,  is  a  function  of  range,  r,  depth,  z,  and  the  azimuthal  angle,  y  .  The  azimuthal 
dependence  of  the  acoustic  field  is  minor  and  ignored  in  the  solution.  Therefore  p  is  a 
function  of  range  and  depth  only,  i.e.  p(r,z). 


Let  the  complex  acoustic  pressure  be  defined  as 


p  1 

(r  ,  z  ,cot) 

1  =  pi 

[r  ,  z  )e  j  a  t 

(3.32) 


and 


(3.33) 
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where  the  /  r  term  accounts  for  azimuthal  spreading  and  u(r,z)  defines  the  two 
/  ylr 

dimensional  acoustic  pressure  field.  Substituting  Eqs.  3.32  and  3.33  into  the  Helmholtz 
equation  yields  the  uncoupled  azimuthal  approximation,  which  is 


(3.34) 


The  following  operators,  denoted  by  the  subscript,  op,  are  defined  as 

(3.35) 


P  =  5 

°p  /a  r 


and 


QoP  -  +  e  +  i 


(3.36) 


where 


£  —  n  —  1 


(3.37) 


and 


(3.38) 


Then  Eq.  3.32  can  be  rewritten  as 


If  u  is  defined  as 


(3.40) 


(3.39) 


then  the  outgoing  wave  can  be  shown  to  satisfy  [Ref.  9] 
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where  the  propagator,  <f>(r) ,  which  is  a  unitary  operator,  steps  the  solution  out  in  range. 
The  split-step  Fourier  method  is  used  to  apply  the  unitary  operator.  The  operator,  Hop,  is 
separated  by  using  the  wide-angle  approximation  [Ref.  10]  of  the  square  root  operation 
such  that 
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Therefore  the  propagator  may  be  expressed  as  [Ref.  1 1] 


The  PE/SSF  method  applies  the  operator 


in  the  kz  domain  where  the  operator,  T  ,  is  defined  as 


II 

1 

\ 

L 

ro 

2 

\ 

l  k°  J 

The  following  convention  is  used  for  the  FFT 


(3.47) 


(3.48) 


(3.49) 


(3.50) 


(3.51) 


(3.52) 
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(3.53) 


v(z)=  FFT  V  (*.-)) 

and 

¥(k!)=  IFFT  V(z)} 


Finally,  the  PE/SSF  method  is  implemented  as 


Ar  /  \ 

/  \  — O  ( z,r+Ar) 

y/(z,r+Ar)=e  2  FFT< 

g  JK^FFz) JJ7J7JI 

~A—Uop(z,r)  .  , 

e  2  W,r)yt 

- 

JJ 

(3.54) 


The  MMPE  model  provides  multiple  outputs  to  the  user.  The  primary  output  for 
use  with  the  OFDM  simulation  is  the  complex  acoustic  pressure.  The  MMPE  model  used 
in  the  simulation  is  a  2D  model  in  range  and  depth  since  the  azimuthal  dependence  was 
neglected  thereby  omitting  the  third  dimension.  The  solution  is  a  grid  of  points  from 
which  the  complex  acoustic  pressure  can  be  evaluated.  The  complex  acoustic  pressure 
output  for  a  given  point  in  the  grid  is  a  function  of  frequency.  The  application  of  this  data 
to  the  OFDM  signal  is  discussed  in  the  receiver  section  covering  the  application  of  the 
acoustic  channel. 

The  source  depth  for  the  shallow  water  (100  m)  acoustic  channels  was  30  m  and 
100  m  for  the  deep  water  (340  m)  acoustic  channels.  The  receiver  depth  varied  with  the 
results  of  the  acoustic  channel  evaluation.  The  receiver  depth  was  chosen  by  determining 
the  depth  with  the  greatest  average  acoustic  pressure  amplitude  across  the  frequency  band 
evaluated.  The  sound  speed  profde  (SSP)  used  for  all  acoustic  channels  is  from  the 
Atlantic  Ocean  just  off  the  eastern  coast  and  is  shown  in  Figure  3.14. 
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D.  RECEIVER 

The  functional  blocks  in  the  receiver  were  presented  in  Fig.  3.1.  For  the  most  part 
the  blocks  of  the  receiver  simply  undo  the  work  of  the  transmitter  to  convert  the 
information  back  to  its  original  form.  Additional  blocks  in  the  receiver  include  the 
channel  estimation  and  equalization  blocks,  additive  white  gaussian  noise  (AWGN), 
synchronization  and  bit  error  analysis. 

1.  AWGN  Addition 

The  simulation  assumes  the  background  noise  present  in  the  channel  is  AWGN. 
The  signal  power  is  calculated  as 
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(3.55) 


X  *2(0 

i  =  1 _ 

K 


where  K  is  the  length  of  the  received  signal.  The  noise  power  to  be  added  is  then 
calculated  as, 


(3.56) 


where  SNR  is  the  desired  SNR  in  dB.  The  received  signal  is  then 


xr\ 

(  > 

n> 

I  =  x\ 

(  > 

n> 

)+w{ 

[n) 

(3.57) 


where  x(n)  is  the  noise  free  signal  at  the  input  to  the  receiver,  which  is  the  output  of  the 
transmitter  and  w(n)  is  AWGN. 


2.  DSB  Demodulator 

The  received  signal  is  first  demodulated  using  a  cosine  signal,  as  presented  in  Eq. 
3.26,  then  the  demodulated  signal  is  filtered  and  decimated.  The  LPF  is  designed  such 
that  spectral  components  above  the  known  carrier  frequency  plus  the  bandwidth  of  the 
transmitted  signal  are  removed.  Then  decimation  by  D  is  performed  to  reduce  the  number 
of  samples  and  the  sampling  rate. 

The  reader  should  note  that  in  the  current  flowchart  of  the  simulation,  some 
processing  of  the  received  signal  occurs  before  the  acoustic  channel  model  is  applied. 
This  is  because  without  source  or  receiver  motion  the  MMPE  model  can  be  applied 
directly  to  the  complex  QAM  symbols  in  each  of  the  OFDM  blocks.  In  reality  the 
channel  effects  would  alter  the  signal  after  transmission  and  before  reception.  In  future 

work  source  and  receiver  motion  will  be  modeled.  At  that  point  the  application  of  the 
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MMPE  model  will  take  place  between  transmission  and  reception.  Therefore  the  blocks 
such  as  the  DSB  Modulator  and  Demodulator  have  been  left  in  the  simulation  even 
though  they  are  not  essential  in  the  current  fonn. 

Signal  detection  theory  is  not  addressed  in  this  thesis.  As  such  it  is  assumed  that 
the  signal  is  detected  by  some  means  and  the  stream  of  data  containing  the  signal  is 
processed.  The  simulation  is  designed  such  that  the  received  signal  of  finite  length 
contains  the  transmitted  signal  of  some  length  less  than  the  length  of  the  received  signal. 
The  information  is  assumed  to  be  received  in  an  AWGN  environment.  The  DSB 
Demodulator  demodulates  the  entire  received  signal  to  baseband.  The  synchronizer  is 
responsible  for  identifying  where  the  information  lies  within  the  received  signal.  Only  the 
portion  of  the  received  signal  containing  information,  i.e.  the  transmitted  signal,  is 
processed  in  the  remainder  of  the  receiver. 


The  demodulating  signal  is  expressed  as 


(3.58) 


where  n  is  defined  from  1  to  the  length  of  the  received  signal,  xr(n).  The  demodulated 
signal  is  then  the  product  of  the  received  signal  with  the  demodulating  signal,  that  is, 


Xr,d  < 

f  \ 

ynj 

l  = 

f  \ 

ynj 

1  xd(n) 

(3.59) 


The  signal  is  baseband  and  real  at  this  point.  Next  the  signal  must  be  decimated  to 
reduce  the  number  of  samples  and  the  sampling  rate  by  the  decimation  factor,  D.  First  the 
signal  must  be  passed  through  an  antialiasing  LPF  to  remove  spectral  components  above 

.  The  LPF  is  again  designed  using  the  Hamming  window  method  with  the  following 
parameters: 
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1.1(2 7tW)  ,  .  1X> 

I  co„  = - =  pass  band  digital  frequency 


F 


s.  mod 


1.25(2^)  ,  .  ir. 

\co,  = - =  stop  band  digital  Irequency 


\Nh  = 


.s.mod 

8  71 
(O  -  CO 


=  filter  order  (number  of  coefficients  is  Nh  + 1 ) 


where  reference  can  be  made  to  Fig.  3.10  for  the  general  filter  characteristics.  The  output 
of  the  filter  is 


Xr,f  1 

(«) 

1  =  h  (n ) 

»*  Xr,J 

[n) 

(3.60) 


where  h(n)  is  the  impulse  response  of  the  LPF  and  *  represents  linear  convolution.  The 
final  step  in  the  DSB  Demodulator  is  decimation  of  the  signal.  Decimation  by  a  factor  D 
results  in  the  expression 


(3.61) 


The  effects  of  the  DSB  Demodulator  are  illustrated  in  Fig.  3.15.  The  signal 
parameters  of  Fig.  3.15  are  the  same  as  those  of  Fig.  3.9.  The  bit  rate  is  5  kbps  using  16 
QAM  and  the  decimation  factor  is  D  is  4.  The  baseband  bandwidth  is  1.501  kHz. 
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Figure  3.15.  DSB  Demodulator  Signal  Processing  Effects. 


3.  Synchronizer 

In  order  to  demodulate  the  subcarriers  an  OFDM  receiver  must  first  perform 
synchronization.  There  are  two  tasks  related  to  synchronization.  The  first  is  determining 
the  location  of  the  OFDM  symbol  blocks.  The  second  is  estimating  and  correcting  for 
frequency  offset  of  the  subcarriers.  In  addition  coherent  receivers  also  require  the  carrier 
phase  of  the  DSB  Modulator  to  be  synchronized.  In  the  simulation  the  phase  of  the  DSB 
Modulator  is  assumed  to  be  synchronized  by  the  OFDM  receiver. 

Recall  from  the  discussion  on  OFDM  theory  that  the  subcarriers  will  remain 
orthogonal  in  the  receiver  only  if  the  transmitter  and  receiver  use  exactly  the  same 
frequencies.  If  not  then  ICI  will  occur.  Practical  oscillators  possess  phase  noise,  which 
results  in  varying  subcarrier  center  frequencies  since  the  frequency  is  the  time  derivative 
of  the  phase.  In  practical  systems  with  non-ideal  oscillators  this  is  an  unavoidable  source 

of  ICI  in  OFDM  systems.  However  Van  Nee  and  Prasad  [Ref.  4]  indicate  that  this  source 
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can  be  minimized  to  acceptable  levels  through  a  well-designed  system.  In  the  simulation 
phase  noise  is  not  modeled  and  therefore  assumed  negligible.  The  sensitivity  of  OFDM 
systems  to  phase  noise  and  frequency  offset  is  one  of  the  major  disadvantages  of  OFDM 
relative  to  single  carrier  systems.  In  single  carrier  systems  phase  noise  and  frequency 
offset  degrade  SNR  but  do  not  produce  interference  [Ref.  4].  Van  Nee  and  Prasad  [Ref. 
4]  provide  a  detailed  discussion  of  these  two  effects  in  OFDM  systems. 

Frequency  offset  is  caused  not  only  by  phase  noise  of  local  oscillators  in  the 
transmitter  but  also  by  Doppler  shift  due  to  non-zero  relative  motion  between  the  source 
and  the  receiver.  The  simulation  assumes  the  transmitter  and  receiver  are  stationary. 
Therefore  frequency  offsets  are  not  addressed  in  this  thesis. 

While  OFDM  systems  are  relatively  sensitive  to  frequency  issues  they  are  quite 
robust  to  timing  errors.  Symbol  timing  may  vary  by  as  much  as  the  duration  of  the  guard 
interval,  tg,  without  producing  ISI  or  ICI.  However  an  optimal  timing  instant  exists  which 
provides  for  maximum  multipath  robustness.  Any  deviation  from  the  optimal  timing 
instant  reduces  the  delay  spread  tolerance  the  system  was  designed  to  handle.  OFDM 
systems  attempt  to  minimize  the  timing  error  relative  to  the  guard  interval. 

The  cyclic  extension  is  used  to  perform  symbol  synchronization  in  the  simulation. 
Recall  that  the  first  tg  samples  received  will  be  the  same  as  the  last  tg  samples  for  each 
OFDM  block  due  to  the  cyclic  extension.  Therefore  correlation  streams  of  data  tg  samples 
long  with  streams  of  the  same  length  but  delayed  by  the  duration  of  the  FFT  interval  will 
identify  the  location  of  the  OFDM  blocks. 

The  synchronization  block  diagram  is  presented  in  [Ref.  4]  and  modified  slightly 
for  the  purposes  of  the  simulation  to  the  fonn  shown  in  Fig.  3.16. 
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Figure  3.16.  Synchronization  Flowchart  [after  Ref.  4]. 


Recall  that  tg  is  the  guard  time  in  continuous  time  and  L  is  the  discrete  time 
number  of  samples  corresponding  to  tg  for  the  applicable  sampling  rate.  Figure  3.17 
illustrates  the  structure  of  an  OFDM  signal  in  serial  form  to  clarify  the  cross  correlation 
used  in  the  synchronizer.  Fig.  3.17  identifies  the  observation  intervals  of  length  N+L 
which  slide  over  the  received  signal  in  discrete  time.  The  first  and  last  L  samples  in  the 
observation  interval  will  be  the  same  when  the  interval  falls  exactly  across  an  OFDM 
symbol  labeled  as  frames  in  Fig.  3. 17.  At  this  point  the  cross  correlation  terms  will  be  the 
greatest.  When  the  observation  interval  includes  samples  from  two  different  frames  such 
that  the  first  and  last  L  samples  are  independent  the  cross  correlation  function  will  have  a 
peak  much  lower  than  when  the  interval  is  exactly  over  one  frame.  In  fact  in  the  limit  as 
the  number  of  samples  over  which  the  cross  correlation  is  performed  is  very  large  the 
ratio  of  the  sidelobe-to-peak  amplitude  will  go  to  zero  [Ref.  4],  Figure  3.18  provides  a 
realization  of  the  correlation  signal  identified  in  Fig.  3.16  for  an  OFDM  signal  with  six 
OFDM  blocks  where  the  first  two  blocks  are  the  channel  estimation  blocks. 


61 


Observation  interval  ODC 


1 

— 

- - - - 1 - 

1  1 

1  1 

• 

i 

i 

i 

• 

Frame  i-1 

Frame  i 

Frame  i+1 

ii 

L 

4  1 

N 

r 

- ► 

s-*- — - — ->v - 

Cyclic  prefix  OFDM  Symbol 

Figure  3.17.  OFDM  Signal  Structure  with  Synchronizer  Observation  Intervals. 


The  peak  finding  algorithm  uses  five  threshold  levels  to  determine  the  existence 
of  peaks  in  the  correlation  signal.  The  threshold  levels  are  adjustable  and  based  on  the 
maximum  value  found  in  the  correlation  signal.  In  the  realization  of  Fig.  3.18  the 
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threshold  levels  are  0.4,  0.45,  0.50,  0.55  and  0.60  of  the  maximum  of  the  correlation 
signal.  The  algorithm  requires  that  for  a  value  to  qualify  as  a  peak  it  must  first  define  an 
inflection  point  in  the  correlation  signal.  Second  the  value  must  register  as  a  peak  at  all 
threshold  levels.  The  impulse  response  and  frequency  response  of  the  acoustic  channel 
model  applied  to  the  transmitted  signal  in  realizing  the  correlation  signal  in  Fig.  3.18  is 
shown  in  Fig.  3.19.  The  impulse  response  of  Fig.  3.19  clearly  indicates  the  existence  of 
three  distinct  receptions  due  to  the  multipath  environment.  Also  note  that  the  frequency 
response  of  Fig.  3.19  illustrates  frequency  selective  fading  and  nearly  linear  phase 
response. 


Acoustic  Channel 


Frequency  Response 


Figure  3.19.  Acoustic  Channel  Filter  Plots  for  Realization  of  Fig.  3.17. 


Van  Nee  and  Prasad  [Ref.  4]  discuss  optimal  timing  by  analyzing  the  correlation 
signal  power.  The  idea  is  to  take  the  samples  that  lie  within  the  OFDM  symbol  that  also 
have  the  greatest  power  and  therefore  the  greatest  SNR.  In  the  simulation  the  signal 
power  is  not  analyzed  since  the  correlator  synchronizer,  while  computationally  intensive, 
was  very  successful  at  identifying  the  boundaries  of  the  OFDM  symbols.  In  most  cases 
the  correct  symbol  boundaries  were  found  to  within  a  few  percent  error  relative  to  the 
length  of  the  cyclic  prefix,  L.  Future  work  might  look  at  the  benefits  of  using  the  optimal 
timing  instant.  Additionally,  the  correlator  synchronizer  used  here  is  computationally 

intensive  and  a  less  demanding  algorithm  would  be  desirable  for  practical  application.  In 
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that  case  the  optimal  timing  instant  would  be  of  greater  value  as  a  less  intense  timing 
synchronization  method  would  likely  be  less  accurate.  A  low  complexity 
synchronization  method  that  seems  to  be  of  value  in  OFDM  is  presented  by  van  de  Beek 
et  al.[Ref.  12].  The  synchronization  is  performed  by  the  Matlab  functions  peakfinder.m 
and  synchronizer.m  contained  in  Appendix  C. 

4.  Composite  Description  of  Serial  to  Parallel  Conversion  through 
FFT/OFDM  Demodulation 

The  serial  to  parallel  (S/P)  converter  takes  the  reduced  length  signal  from  the 
synchronizer  and  converts  it  to  parallel  fonn.  The  resulting  matrix  of  samples  will  have 
the  number  of  rows  corresponding  to  the  block  length  determined  by  the  synchronizer 
and  the  number  of  columns  will  be  the  number  of  OFDM  symbols  found  by  the 
synchronizer.  Recall  that  each  peak  in  the  correlation  signal  corresponds  to  the  start  of  an 
OFDM  symbol.  The  matrix  of  real  time  samples  from  the  S/P  converter  contain  the  cyclic 
prefix  which  is  removed  from  each  OFDM  symbol,  i.e.  column  in  the  matrix,  by  the 
cyclic  extension  removal  block. 

The  OFDM  Demodulator  simply  performs  an  N  point  FFT  of  the  matrix  from  the 
cyclic  extension  removal  block.  The  FFT  can  be  expressed  as 

(3.62) 


(3.63) 

where  W_N  e  C  VxA  and  xN  e  RNxM  therefore  XN  e  CVxV/ .  The  resulting  matrix  of 
complex  samples  represent  the  transmitted  complex  QAM  symbols.  The  OFDM  symbols 
making  up  each  column  of  the  matrix  are  conjugate  symmetric  and  only  the  active 
subcarriers  are  of  interest  for  receiving  the  information  in  the  signal. 


The  FFT  can  also  be  expressed  in  matrix  form  as 
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5. 


Acoustic  Channel  Application 


The  output  of  the  MMPE  model  used  in  the  simulation  is  in  the  form  of  a  complex 
transfer  function.  The  transfer  function  is  applied  directly  to  the  complex  QAM  symbols 
present  after  the  OFDM  demodulator.  The  generic  output  from  the  MMPE  model  is  a 
grid  of  points  in  depth  and  range.  The  desired  output  is  chosen  for  a  given  point.  The 
complex  transfer  function  represents  the  complex  acoustic  pressure  as  a  function  of 
frequency  present  at  the  chosen  point  in  the  grid. 

The  transfer  function  contains  values  for  the  frequencies  between  fc  and  fc  + 

NA fsc  .  This  band  represents  the  frequencies  that  would  lie  in  the  real  frequency  region 

of  the  baseband  signal  corresponding  to  the  signal  modulated  to  fc.  Recall  that  only  the 
active  subcarriers  carry  non-zero  signal  power  and  therefore  determine  the  bandwidth  of 
the  transmitted  signal.  In  order  to  apply  the  transfer  function  to  the  complex  QAM 
symbols  it  must  be  made  conjugate  symmetric  by  application  of  Eq.  3.12.  The  complex 
QAM  symbols  after  application  of  the  transfer  function  is 


Y  =  X  .*  H 


(3.64) 


where  le  CNx iNbhet+0)  and  H_e  C  Vx  therefore  7e  cNx(Nbhct+2) .  The  variable  X 
represents  the  matrix  of  complex  QAM  symbols  after  the  OFDM  demodulator,  H 

represents  the  complex  transfer  function  and  Y.  is  the  output.  The  .*  in  Eq.  3.64 

V 

indicates  array  multiplication  such  that  each  column  of  _  is  weighted  by  the  vector 

H  .  The  number  of  columns  in  XandY  is  N hlock  +  2  due  to  the  2  blocks  of  complex 
QAM  symbols  used  for  channel  estimation  and  equalization. 


6.  Channel  Equalization 

The  channel  equalization  scheme  simply  attempts  to  invert  the  effects  of  the 

channel  such  that  the  complex  QAM  symbols  after  equalization  will  be  the  same  as  those 
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that  were  generated  in  the  transmitter.  The  first  two  OFDM  symbols  are  used  to  estimate 
the  channel  transfer  function.  The  generation  of  the  channel  estimation  blocks  was 
discussed  in  the  QAM  Modulator  section.  The  channel  estimate  is 


(3.65) 

where  Xr  represents  the  received  complex  QAM  symbols  in  the  channel  estimation 
blocks,  X  is  the  complex  QAM  symbols  present  in  the  QAM  Modulator  of  the  transmitter 
and  Hest  is  an  A  by  2  matrix  of  complex  values  representing  estimates  of  the  magnitude 
and  phase  of  the  channel  transfer  function  at  frequencies  corresponding  to  the  subcarrier 
center  frequencies.  Array  division  is  implied  in  Eq.  3.65.  The  values  of X are  assumed  to 
be  known  at  the  receiver.  Each  row  of  Hest  is  averaged  to  give  a  vector  of  values 
representing  the  estimate  of  the  channel  transfer  function.  That  is, 

(3.66) 

In  the  simulation  the  estimates  of  the  channel  transfer  function  are  performed  over 
the  active  subcarriers  only  rather  than  over  all  N  as  using  the  zero  subcarriers  would 
result  in  division  by  zero  in  Hest. 

Finally  the  complex  QAM  symbols  are  equalized  by  defining 

(3.67) 

where  the  division  is  array  division  such  that  the  columns  of  Xr  are  weighted  by  the 
inverse  of  Hest . 

A  realization  of  the  channel  estimates  is  shown  in  Fig.  3.20.  The  magnitude  of  the 
channel  estimate  for  the  two  channel  estimation  blocks  and  the  mean  of  the  two  estimates 
is  plotted  versus  the  FFT  bin  number  which  is  analogous  with  subcarrier  number  or 
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frequency.  The  plot  is  discontinuous  at  high  and  low  bin  numbers  due  to  division  by  zero 
that  would  result  where  the  zero  symbols  are  interspersed  to  evenly  shape  the  power 
distribution  as  discussed  in  the  QAM  Modulator  section.  The  plots  are  over  the  active 
subcarriers  only.  Note  that  Hmean  in  Fig.  3.20  refers  to  the  mean  of  the  two  channel 
estimates  Hsys(l)  and  Hsys(2). 


Figure  3.20.  Channel  Estimates,  7V=2048,  16  QAM,  5  kbps,/c  =  6  kHz,  N  act  =  750  . 


7.  QAM  Demodulator 

The  QAM  demodulator  transforms  the  received  complex  QAM  symbols  to  points 
on  a  Q-av y  QAM  constellation.  Therefore  the  output  will  be  integers  from  0  to  Q-l. 

The  first  step  is  to  strip  the  first  N act  +  1  symbols  from  each  block  as  these  are 
the  symbols  that  carry  the  transmitted  information.  This  implies  that  the  receiver  knows 
the  number  of  active  subcarriers.  This  information  is  assumed  to  be  contained  in  protocol 
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information  that  would  be  carried  as  overhead  in  the  transmission  system.  Recall  that  the 
zero  symbols  that  were  added  in  the  QAM  modulator  are  for  data  management  purposes 
and  could  just  as  well  be  protocol  information.  Protocol  is  not  addressed  in  this  thesis  but 
is  assumed  to  be  present  in  making  some  assumptions  that  are  required  in  the  receiver. 

The  zero  symbols  are  removed  and  the  magnitude  and  phase  of  each  symbol  is 
determined  and  converted  to  Cartesian  coordinates.  The  real  and  imaginary  components 
of  the  received  symbols  correspond  to  the  in-phase  and  Quadrature  components 
respectively.  The  operation  of  the  QAM  Demodulator  up  to  this  point  is  illustrated  in  Fig. 
3.21. 


Figure  3 .2 1 .  QAM  Demodulator. 


The  next  step  is  to  map  the  in-phase  and  Quadrature  components  to  Q-ary  QAM 

constellation  points.  The  mapping  of  constellation  points  to  in-phase  and  Quadrature 

components  was  explained  in  the  QAM  modulator.  In  that  case  the  constellation  point 
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was  exact  and  resulted  in  integer  values  for  the  in-phase  and  Quadrature  components. 
However  in  mapping  back  to  QAM  constellation  points  the  in-phase  and  Quadrature 
components  will  not  in  general  be  integer  values.  Therefore  an  algorithm  is  required  to 
determine  the  constellation  point  to  which  the  received  symbol  most  closely  corresponds. 

This  operation  is  performed  by  the  Matlab  function  qaskdecomod.m  contained  in 
Appendix  C.  A  2  by  Q  matrix  containing  the  in-phase  and  Quadrature  values 
corresponding  to  each  point  in  a  £>-ary  QAM  constellation  is  created.  The  in-phase  and 
Quadrature  components  of  the  received  symbols  are  compared  to  the  vectors  in  the 
constellation  matrix  to  determine  the  constellation  point  that  is  nearest  to  the  received 
symbol.  The  received  symbol  is  mapped  to  this  constellation  point  by  its  corresponding 
integer  value.  The  process  is  illustrated  in  Fig.  3.22.  The  example  in  Fig.  3.22  is  for  a 
constellation  size  of  4  for  simplicity  although  4-ary  QAM  constellations  are  not  used  for 
data  transmission  in  the  simulation.  The  error  analysis  is  performed  for  each  received 
symbol  resulting  in  a  mapping  of  an  integer  from  0  to  Q- 1  for  each  in-phase  and 
Quadrature  component  pair.  Therefore  the  output  of  the  QAM  Demodulator  is  a  vector  of 

integer  values  of  length  N  act  *  N bi0ck  N0 . 
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8.  Composite  Description  of  Block  De-Interleaving  through  Error 
Analysis 

The  Block  De-interleaver  simply  undoes  the  interleaving  of  data  performed  by  the 
Block  Interleaver  in  the  transmitter.  It  is  assumed  that  the  receiver  knows  the  code  word 
length  of  the  Reed  Solomon  code  as  the  code  vector  supplied  to  the  deinterleaver  is 
reshaped  into  a  matrix  with  n  columns  in  preparation  for  Reed  Solomon  decoding. 

The  Reed  Solomon  decoder  applies  an  algorithm  in  the  Matlab  communications 
toolbox,  rsdeco.m,  in  decoding.  Recall  the  code  word  length  is  n  and  the  message  length 
is  k.  Therefore  the  input  matrix  has  n  columns  and  the  output  matrix  has  k  columns  where 
each  row  corresponds  to  a  codeword  when  coded  and  a  message  when  decoded. 

The  integers  representing  QAM  symbols  must  be  converted  to  binary  to  continue 
the  signal  reconstruction.  Each  integer  is  converted  to  a  binary  word  with  q  bits  where  q 
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is  the  number  of  bits  per  QAM  symbol.  The  output  of  this  block  is  a  vector  of  binary 
digits. 

The  quantization  buffer  simply  reshapes  the  bit  stream  to  a  matrix  with  eight 
columns,  since  8-bit  quantization  was  performed  in  the  transmitter.  As  stated  earlier  the 
level  of  quantization  is  adjustable  and  detennined  by  the  variable  bits  in  the  transmitter. 
In  general  the  bit  stream  is  reshaped  to  a  bits  column  matrix  in  preparation  for  signal 
reconstruction. 

Signal  reconstruction  consists  of  two  steps.  The  first  is  the  conversion  of  the 
binary  words  of  length  bits  to  integers.  Second  the  integer  values  are  mapped  to  the 
corresponding  value  of  the  information  signal  in  the  transmitter.  Refer  back  to  Fig.  3.4 
which  is  the  illustration  of  the  quantization  process. 

The  BER  is  determined  by  comparing  the  transmitted  and  received  bit  streams. 
An  error  occurs  when  the  received  bit  is  different  than  the  corresponding  transmitted  bit. 
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IV.  RESULTS 


The  results  of  the  thesis  are  divided  into  2  major  sections.  The  first  section 
develops  the  theoretical  performance  of  an  OFDM  system.  The  second  presents  the 
results  of  the  OFDM  simulation  for  several  combinations  of  system  and  acoustic  channels 
parameters. 

A.  SYSTEM  THEORETICAL  PERFORMANCE 

The  theoretical  perfonnance  of  the  system  is  first  analyzed  by  considering  a  single 
carrier  system.  An  OFDM  system  is  equivalent  to  a  set  of  independent  and  ISI  free  single 
carrier  QAM  systems.  Therefore  a  large  part  of  the  analysis  of  an  OFDM  system  can  be 
performed  by  analyzing  a  single  carrier  QAM  system.  The  single  carrier  analysis  is 
extended  to  the  multicarrier,  OFDM,  system  to  complete  the  analysis. 


1.  Single  Carrier  Analysis 


by 


The  energy  per  two-dimensional  symbol  in  a  square  QAM  constellation  is  given 


(4.1) 


where  d  is  the  distance  between  points  in  the  QAM  constellation  and  Q  is  the  size  of  the 
QAM  constellation.  All  points  in  the  constellation  are  equally  likely  and  the  constellation 
is  centered  at  the  origin  and  therefore  has  zero  mean.  Equation  4.1  is  exact  only  for 
square  QAM  constellations  but  is  still  a  good  approximation  for  the  average  transmit 
energy  for  non-square  QAM  constellations  [Ref.  1].  For  an  ISI  free  channel  with  gain, 
\h\  ,  the  probability  of  two  dimensional  symbol  error  is  approximately 


(4.2) 
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where  dmm  is  the  minimum  distance  between  QAM  constellation  points  at  the  output  of 
the  channel  and  Qerj(x)  is  the  error  function,  each  of  which  are  defined  as  [Ref.  3] 


and 


(4.3) 


(4.4) 


The  probability  of  symbol  error  per  dimension,  which  is  PJ 2,  for  communication 
systems  should  be  less  than  10‘4.  Therefore  the  following  requirement  must  be  met, 


U . 


v  2  o  j 


11.8  {dB)+r,„{dB)-rc(dB) 


(4.5) 


The  variable  ym  is  the  margin,  which  provides  for  additional  performance  in  order 
to  overcome  unpredictable  channel  degradations  and  yc  is  the  coding  gain.  Note  that 
when  ym  =  yc  =  0  the  system  is  uncoded  and  has  no  margin.  Therefore  11.8  dB  is  the 

argument  required  in  the  Q(.)  function  to  get  =  10  4 .  Coding  increases  dmm  and 

therefore  reduces  the  1 1.8  dB  by  the  amount  of  coding  gain.  Using  a  margin  increases  the 
1 1.8  dB  required  by  the  amount  of  margin  desired.  [Ref.  3] 

Equation  4.1  can  be  rewritten  as 


i 

b 

ro 

H 

2 

KS  —  i  1  ~ 

d2 

min 

(4.6) 


which  leads  to  the  definition  of  SNR  gap,  or  normalized  SNR,  T , 
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(4.7) 

where  o1  is  the  variance  or  power  of  the  AWGN. 

Therefore,  for  the  target  of  =  10  4 ,  the  following  expression  can  be  written 
for  the  SNR  gap, 

7>l{dB)  =  U.%+Ya-Yc{dB)\  (48) 

Then  by  taking  the  base  2  log  of  Eq.  4.6  and  substituting  into  Equation  4.7  for 

dmin,  the  number  of  bits  that  can  be  carried  by  QAM  at  the  target  PJ2  is  calculated  as 

I  ,  (r>\  .  (,  SNR) 

bp,  =  log  2{Q)  =  log2l  1  +  —  I  (49) 

where  the  channel  output  SNR  is  defined  as 

(4.10) 

Note  that  Eq.  4.9  is  the  Shannon-Hartley  theoretical  channel  capacity  with  the 
SNR  reduced  by  a  factor  of  the  SNR  gap.  As  the  SNR  gap  goes  to  1  (0  dB)  the  achievable 
data  rate  of  the  system  approaches  capacity.  As  can  be  seen,  the  SNR  gap  is  a  measure  of 
the  degradation  of  the  system  from  optimal.  The  number  of  bits  that  can  be  carried  by  a 
QAM  system  as  computed  using  Eq.  4.9  will  not  be  an  integer  but  is  typically  rounded 
down  to  the  nearest  integer.  If  a  square  QAM  constellation  is  desired  then  it  is  rounded 
down  to  the  nearest  even  integer.  The  SNR  gap  for  any  QAM  system  with  probability  of 
one  dimensional  symbol  error  equal  to  10'4  is 

r  =  7.029  +  ym  -  y c  r4.1i) 
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2. 


OFDM  (Multi-Carrier)  Analysis 


With  the  single  carrier  analysis  complete,  multi-carrier  analysis  can  now  be 
performed  as  an  extension  of  the  single  carrier  results.  In  a  generic  multicarrier  system 
the  probability  of  error  is  an  average  of  the  probability  of  error  for  each  of  the 
subchannels.  Therefore  the  overall  probability  of  error  will  be  dominated  by  the  weakest 
subcarriers.  [Ref.  3] 


This  weakness  was  pointed  out  earlier  in  the  thesis.  Forward  error  correction 
coding  with  block  interleaving  is  used  to  average  out  the  probability  of  error  over  the 
active  subcarriers.  Therefore  in  the  following  analysis  the  probability  of  error  is  assumed 
equal  over  all  subchannels.  Directly  applying  the  single  carrier  analysis  yields  the 
following  results  for  the  ith  subchannel 


>4 

1 

1 

3 

p 

*** . 

1 

1 

H  , 

2  J  2 

d  i 

1  0  — 
4(7  r 

4(7  ,2 

(4.12) 


where  the  subscript  i  denotes  any  values  that  are  subchannel  dependent.  The  maximum 
number  of  bits  that  can  be  carried  on  the  i'h  subchannel  with  a  given  margin  and  coding 
gain  is 


(4.13) 


where 


(4.14) 


The  subsymbol  energy,  ,  is  held  constant  over  all  subchannels  that  are  active 

and  is  zero  on  the  inactive  subchannels.  The  water  pouring  distribution  (see  Proakis  [Ref. 

13]  or  Gallager  [Ref.  14]),  is  a  better  energy  distribution  but  Cioffi  has  determined  that 
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the  on/off  energy  distribution  is  very  close  to  optimal  yet  is  much  easier  to  compute  and 
implement.  [Ref.  3] 

Figure  4.1  plots  curves  of  bpe  versus  received  channel  SNR  for  three  different 
values  of  one-dimensional  probability  of  symbol  error,  Pej.  Figure  4.1  assumes  that  the 
coding  gain  is  equal  to  the  margin  which  may  be  zero.  Therefore  actual  results  in 
simulation  or  that  may  be  seen  in  real  world  application  may  shift  left  or  right  as  a 
function  of  the  margin  and  coding  gain.  If  the  coding  gain  is  greater  than  the  margin  the 
curves  will  shift  left  thereby  predicting  greater  numbers  of  bits  per  subcarrier  at  a  given 
received  SNR.  The  coding  gain  achieved  in  the  simulation  is  not  calculated  but  is 
assumed  to  be  around  6  to  8  dB  which  will  be  shown  to  be  accurate  in  the  following 
results. 


Figure  4. 1 .  Channel  Bit  per  Subcarrier  Capacity. 


The  bit  rate,  R,  of  the  total  system  is  then  simply  the  product  of  the  active 
subcarriers  with  the  number  of  bits  per  subcarrier  divided  by  the  OFDM  symbol  period, 
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T.  Since  all  active  subcarriers  have  the  same  number  of  bits,  bsc,  which  is  rounded  down 
to  the  nearest  integer  below  bp  A  ,  the  system  bit  rate  is 


Figure  4.2  plots  system  bandwidth  versus  bit  rate  for  several  sizes  of  QAM 
constellations. 


Figure  4.2.  Bandwidth  vs  Bit  Rate  for  Q-ary  QAM  OFDM  Signals. 


B.  SIMULATION  RESULTS 

The  remainder  of  this  chapter  presents  the  results  of  the  simulation  over  the 
acoustic  channel  models  developed.  The  main  parameters  of  the  acoustic  channel  are  the 
water  depth,  the  bottom  roughness  in  meters  of  standard  deviation,  the  range  between  the 
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receiver  and  the  transmitter  and  the  frequency  band  over  which  the  model  is  evaluated. 
The  water  depths  evaluated  are  either  100  m  or  340  m.  Three  different  bottom  roughness 
values  of  0,  2  and  4  m  standard  deviation  were  evaluated.  The  ranges  evaluated  were 
either  2  km  or  4  km. 

Two  of  the  available  forms  of  output  from  the  MMPE  model  are  included  for 
every  acoustic  model  evaluated.  The  first  is  a  plot  of  transmission  loss  as  a  function  of 
frequency  and  depth.  The  second  is  a  plot  of  transmission  loss  versus  depth  and  time. 
These  plots  show  the  ray  paths  present  in  the  water  column  over  time  such  that  for  a 
given  depth  the  multipath  arrivals  are  apparent. 

While  each  acoustic  model  was  evaluated  for  ranges  of  2  and  4  km,  only  the  plots 
for  2  km  ranges  are  provided.  Figure  4.3  shows  the  transmission  loss  for  the  acoustic 
channel  that  is  100  m  deep  and  has  a  bottom  roughness  of  4  m.  The  effect  of  the  bottom 
roughness  on  the  acoustic  energy  is  clearly  seen  in  Figure  4.3  and  can  be  compared  with 
the  transmission  loss  curves  of  the  same  acoustic  channel  with  less  bottom  roughness, 
which  are  Figures  B.l  and  B.3,  to  further  illustrate  the  effect.  Note  the  interference 
patterns  and  blending  of  transmission  near  the  interface  that  is  much  less  pronounced  as 
the  bottom  roughness  goes  down.  Figure  4.4  shows  the  acoustic  energy  arrival  as  a 
function  of  time  and  depth  for  the  same  acoustic  channel.  The  severe  multipath 
environment  present  in  shallow  water  acoustic  channels  is  readily  apparent  in  Figure  4.4. 
All  other  acoustic  channel  plots  are  found  in  Appendix  B. 

The  simulation  was  performed  for  several  combinations  of  bit  rate,  center 
frequency  and  constellation  sizes.  Each  combination  was  simulated  10  times  and  then  the 
mean  of  the  result  was  taken.  The  primary  output  of  the  simulation  and  this  thesis  are 
plots  of  BER  versus  received  SNR.  Recall  that  the  simulation  is  performed  with  the 

margin,  7m  ,  set  equal  to  zero.  Accepted  values  for  7m  range  from  6  to  12  dB,  with 
Cioffi  [Ref.  3]  recommending  6  dB. 
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Transmission  Loss  (dB  re  lm) 
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Figure  4.3.  Transmission  Loss,  Bottom  Roughness  4  m,  Water  Depth  100  m. 

The  margin  can  be  considered  as  an  adjustment  for  simulations  to  shift  the  results 

towards  what  can  be  expected  in  reality.  Recall  also  that  the  coding  gain,  Yc  ,  is  estimated 

to  be  4  to  6  dB.  Therefore,  since  the  coding  gain  is  included  in  the  results  of  the 
simulation,  in  order  to  check  the  results  versus  the  theoretical  prediction  of  Fig.  4.1  the 
BER  versus  received  SNR  must  be  shifted  to  the  right  by  approximately  1 1  dB. 
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Figure  4.4.  Acoustic  Energy  Arrival,  Bottom  Roughness  4  m,  Water  Depth  100  m. 


The  raw  results  are  provided  in  Appendix  A.  In  these  plots  each  data  point 
indicates  the  result  of  one  run  of  the  simulation.  The  solid  red  line  is  the  mean  of  the  10 
runs  performed.  The  composite  results  are  provided  in  Figures  4.5  through  4.8  for  QAM 
constellation  sizes  of  8,  16,  32  and  64.  Each  solid  line  in  the  figures  is  the  mean  of  the 
corresponding  series  of  10  runs  plotted  individually  in  Appendix  A.  Figure  4.5  presents 
the  results  for  all  evaluations  using  8-ary  QAM. 
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8  QAM  Composite  Results 


Figure  4.5.  8  QAM  Composite  Results. 


Figure  4.6.  16  QAM  Composite  Results. 
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Figure  4.7.  32  QAM  Composite  Results. 


Figure  4.8.  64  QAM  Composite  Results. 
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Figure  4.9  provides  an  appreciation  for  the  nature  of  the  decoded  signal  and  its 
frequency  spectrum  with  respect  to  the  information  signal  and  its  frequency  spectrum. 


Transmitted  Information  Signal 


Figure  4.9.  16  QAM  Signal  Realization,/;  =  6  kHz,  R  =  5  kbps,  SNR  =  5.3  dB 

BER  =  1.6  E-2. 


The  BER  is  unacceptable  by  communication  system  standards.  The  information 
signal  for  the  simulation  is  a  voice  stream.  Playback  of  the  received  signal  at  the  BER  of 
Figure  4.9  contains  minor  static  but  is  clearly  understandable.  The  zoomed-in  subplots 
show  the  deviation  of  the  received  information  signal  spectrum  from  the  transmitted 
information  signal  spectrum.  One  can  see  that  where  there  is  significant  power  in  the 
frequency  spectrum  of  the  transmitted  information  signal  the  received  spectrum  is  very 
near  that  of  the  transmitted  spectrum.  Comparing  the  time  domain  plots  of  the  transmitted 
and  received  information  signals  one  can  see  the  deviations  of  the  received  signal  relative 
to  the  original  information  signal.  Figure  4.10  shows  the  received  QAM  constellation  for 
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the  realization  of  Figure  4.9  after  channel  estimation  and  equalization  but  prior  to  Reed 
Solomon  decoding. 


Received  QAM  symbols  prior  to  Reed  Solomon  decoding 
16  QAM,  BER=  1.6  E-2 
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Figure  4.10.  16  QAM  Signal  Realization,^  =  6  kHz,  R  =  5  kbps,  SNR  =  5.3  dB 


BER  =  1.6  E-2. 


Figures  4.11  and  4.12  provide  a  contrast  to  the  realizations  of  Figures  4.9  and  4.10 
at  a  greater  received  SNR  of  7.4  dB  and  the  resulting  smaller  BER  of  2.5  E-4.  The  BER 
is  still  greater  than  the  acceptable  limit  of  1  E-4.  However  one  can  see  that  the  received 
signal  is  essentially  an  exact  replica  of  the  transmitted  information  signal.  Also  the 
powerful  effect  of  the  FEC  can  be  seen  from  the  noisy  constellations  of  Figures  4.10  and 
4.12. 
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BER  =  2.5  E-4. 


Note  that  the  zoomed  plots  of  the  frequency  spectrums  of  Fig.  4.11  show  that  the 
received  spectra  is  essentially  the  same  as  that  of  the  transmitted  information  signal 
spectrum  even  in  the  areas  of  the  spectrum  where  the  transmitted  information  signal  has 
relatively  low  power. 
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Received  QAM  symbols  prior  to  RS  decoding 
16  QAM,  BER  2.5  E-4 
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Figure  4.12.  16  QAM  Signal  Realization,/,  =  6  kHz,  R  =  5  kbps,  SNR  =  7.4  dB 

BER  =  2.5  E-4. 
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V.  CONCLUSIONS 


The  primary  thrust  of  this  thesis  was  to  develop  a  computer  simulation  of  an 
OFDM  system  for  underwater  acoustic  communication  using  a  PE  based  model  of  the 
acoustic  channel.  The  simulation  was  validated  by  comparison  of  the  results  of  the 
simulations  with  the  known  theoretical  results  established  in  this  thesis  and  those  from 
ocean  experiments.  The  comparison  of  the  results  with  what  is  expected  from  the  theory 
was  made  in  conjunction  with  presenting  the  results.  Coatelan  and  Glavieux  [Ref.  5] 
experimented  with  a  OFDM  communication  system  in  shallow  water  acoustic  channels 
with  results  as  shown  in  Table  5.1.  The  QAM  constellation  size  used  by  Coatelan  and 
Glavieux  [Ref.  1]  is  not  known.  Note  that  the  water  depth  of  the  experiments  is  much 
shallower  than  for  the  acoustic  channels  used  in  the  simulations  of  this  thesis.  The  FEC 
code  employed  in  the  experiments,  which  uses  rate  A  convolutional  coding  with  a 
constraint  length  of  seven  as  well  as  a  soft  Viterbi  decoder  working  on  a  64  state  trellis 
[Ref.  1],  is  more  advanced  than  that  of  the  simulation,  however  the  coding  gain  is  likely 
within  a  couple  of  dB  of  the  coding  gain  for  the  simulation  as  more  than  a  few  dB  of 
coding  gain  improvement  beyond  5  dB  is  difficult  to  obtain.  Comparing  the  composite 
results  presented  earlier  for  the  8  and  16  QAM  constellations  and  applying  a  margin  of  6 
to  12  dB  one  can  see  that  simulation  results  are  comparable  to  the  experimental  results  of 
Table  5.1. 


Transmission 

distance 

Water 

depth 

Seafloor 

type 

Emitted 

number 
of  bit 

BER  after 
demodulation 

BER  after 
decoding 

estimated 

average 
level  Et,/N0 

650  m 

25  m 

sludge 

18500 

3,7. 10"3 

<  1.10"4 

27,5  dB 

1920  m 

27  m 

sand 

18500 

2,4.1  O'2 

<  1.1  O'4 

16  dB 

Table  1.  Multicarrier  transmission  performances  during  two  stations 

Table  5.1.  Experimental  Multicarrier  Results  [from  Ref.  1]. 
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The  benefits  of  the  work  are  the  verification  of  the  feasibility  of  using  OFDM  as  a 
method  of  underwater  acoustic  communication  and  to  provide  a  test-bed  for  future  work 
to  develop  improved  methods  for  the  blocks  that  make  up  the  OFDM  system. 

There  is  significant  opportunity  for  future  work  in  this  area.  Many  of  the 
algorithms  applied  in  the  OFDM  system  are  simple  and  could  be  improved  to  enhance  the 
performance  of  the  OFDM  system. 

The  modulation  is  DSBSC  which  uses  twice  the  bandwidth  of  an  Upper  Side 
Band  (USB)  system.  There  are  no  obvious  reasons  why  USB  modulation  could  not  be 
used  to  improve  the  spectral  efficiency.  This  idea  also  has  benefit  in  working  on  a  multi¬ 
user  application  of  OFDM. 

Communication  system  protocols  must  be  established  to  make  the  system  useful 
for  real  world  application  and  appears  to  have  been  neglected  so  far  in  the  area  of 
underwater  acoustic  communication.  The  IEEE  standard  802.11  and  the  European  digital 
terrestrial  TV  broadcast  uses  OFDM  and  has  established  protocols  that  may  be  applicable 
to  underwater  communication  systems. 

Synchronization  algorithms  have  significant  merit  for  future  study  and 
improvement.  This  study  does  not  involve  Doppler  effects  since  it  assumes  that  the 
transmitter  and  receiver  are  stationary. 

The  application  of  the  channel  transfer  function  must  be  modified  in  modeling 
source  or  receiver  motion.  Frequency  offset  due  to  wave  action  and  phase  noise  are  not 
simulated  either.  Both  of  these  are  real  concerns  in  development  of  an  actual  system  and 
are  worthy  of  simulation  and  analysis.  The  synchronization  method  must  be  modified  in 
the  simulation  to  account  for  phase  issues.  In  addition  the  synchronization  method  is  very 
computationally  intensive  and  simplification  of  the  system  would  benefit  the  simulation 
and  any  resulting  real  system. 

Channel  estimation  and  equalization  assumes  a  linear  time  invariant  acoustic 
channel.  The  channel  is  known  to  be  time  varying  with  the  rate  of  variance  usually  low 
but  still  of  concern.  Therefore  a  method  of  updating  the  channel  estimation  is  necessary. 
Also  while  the  estimation  method  used  in  the  simulation  is  very  robust  there  may  be 
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better  methods  that  use  less  bandwidth.  FEC  coding  applied  is  elementary  and  can  be 
improved. 

Thus  far  work  to  improve  the  system  and  the  degree  of  reality  of  the  simulation 
have  been  discussed.  Of  equal  importance  is  the  study  of  the  key  parameters  on  the 
system  perfonnance.  The  peak  to  average  power  ratio  of  the  transmit  signal  is  of  concern 
in  OFDM  and  has  not  been  addressed  in  this  thesis. 
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APPENDIX  A.  SIMULATION  RESULTS 


The  figures  in  this  appendix  represent  the  results  of  the  simulation  for  all  of  the 
acoustic  channels  evaluated.  The  parameters  of  the  system  are  provided  in  each  figure  as 
well  as  the  main  channel  parameters.  Each  figure  includes  a  zoomed  inset  which 
illustrates  the  BER  performance  near  the  target  level  of  10'4. 


Figure  A.  1 .  BER  vs  SNR  for  8  QAM,  fc  =  8  kHz,  R  =  5  kbps,  Range  =  2  km,  Rough 

=  Om,  Water  Depth  =  100m. 
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Figure  A.2.  BER  vs  SNR  for  8  QAM,  fc  =  8  kHz,  R  =  5  kbps,  Range  =  4  km,  Rough 

=  Om,  Water  Depth  =  100m. 


SNR  (dB) 


Figure  A. 3.  BER  vs  SNR  for  8  QAM,  fc  =  10  kHz,  R  =  6  kbps,  Range  =  4  km,  Rough 

=  0m,  Water  Depth  =  100m. 
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Figure  A. 4.  BER  vs  SNR  for  8  QAM,  fc  =  10  kHz,  R  =  6  kbps,  Range  =  2  km,  Rough 

=  0m,  Water  Depth  =  100m. 


Figure  A. 5.  BER  vs  SNR  for  8  QAM,  fc  =  8  kHz,  R  =  5  kbps,  Range  =  2  km,  Rough 

=  0m,  Water  Depth  =  100m. 
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Figure  A. 6.  BER  vs  SNR  for  8  QAM,  fc  =  12  kHz,  R  =  7.5  kbps,  Range  =  4  km, 

Rough  =  Om,  Water  Depth  =  100m. 


Figure  A. 7.  BER  vs  SNR  for  16  QAM,  fc  =  6  kHz,  R  =  5  kbps,  Range  =  2  km,  Rough 

=  0m,  Water  Depth  =  100m. 


96 


Figure  A. 8.  BER  vs  SNR  for  16  QAM,  fc  =  6  kHz,  R  =  5  kbps,  Range  =  2  km,  Rough 

=  2m,  Water  Depth  =  100m. 


Figure  A.. 9  BER  vs  SNR  for  16  QAM,  fc  =  6  kHz,  R  =  5  kbps,  Range  =  2  km,  Rough 

=  4m,  Water  Depth  =  100m. 
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Figure  A.  10.  BER  vs  SNR  for  16  QAM,  fc  =  6  kHz,  R  =  5  kbps,  Range  =  4  km,  Rough 

=  0m,  Water  Depth  =  100m. 


SNR  =  21dB,  f  =  6kHz,  R  =  5kbps,  Range  =  4  km.  Rough  =  4,  Water  Depth  =  100  m 


SNR  (dB) 


Figure  A.  1 1 .  BER  vs  SNR  for  16  QAM,  fc  =  6  kHz,  R  =  5  kbps,  Range  =  4  km,  Rough 

=  4m,  Water  Depth  =  100m. 
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Figure  A.  12.  BER  vs  SNR  for  16  QAM,  fc  =  6  kHz,  R  =  5  kbps,  Range  =  2  km,  Rough 

=  Om,  Water  Depth  =  340m. 


Figure  A.  13.  BER  vs  SNR  for  16  QAM,  fc  =  6  kHz,  R  =  5  kbps,  Range  =  4  km,  Rough 

=  2m,  Water  Depth  =  100m. 
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Figure  A.  14.  BER  vs  SNR  for  16  QAM,  fc  =  6  kHz,  R  =  5  kbps,  Range  =  2  km,  Rough 

=  2m,  Water  Depth  =  340m. 


SNR  =  21dB,  f  =  6kHz,  R  =  5kbps,  Range  =  4  km,  Rough  =  0,  Water  Depth  =  340  m 


SNR  (dB) 


Figure  A.  15.  BER  vs  SNR  for  16  QAM,  fc  =  6  kHz,  R  =  5  kbps,  Range  =  4  km,  Rough 

=  0m,  Water  Depth  =  340m. 
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Figure  A.  16.  BER  vs  SNR  for  16  QAM,  fc  =  6  kHz,  R  =  5  kbps,  Range  =  4  km,  Rough 

=  2m,  Water  Depth  =  340m. 


Figure  A.  17.  BER  vs  SNR  for  16  QAM,  fc  =  12  kHz,  R  =  10  kbps,  Range  =  2  km, 

Rough  =  0m,  Water  Depth  =  100m. 
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SNR  =  2 1  dB,  f  =  1 2kHz,  R  =  1  Okbps,  Range  =  4  km.  Rough  =  0,  Water  Depth  =  1 00  m 


SNR  (dB) 


Figure  A.  18.  BER  vs  SNR  for  16  QAM,  fc  =  12  kHz,  R  =  10  kbps,  Range  =  4  km, 

Rough  =  0m,  Water  Depth  =  100m. 


Figure  A.  19.  BER  vs  SNR  for  16  QAM,  fc  =  8  kHz,  R  =  6.67  kbps,  Range  =  2  km, 

Rough  =  0m,  Water  Depth  =  100m. 
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SNR  =  21dB,  f  =  8kHz,  R  =  6.6667kbps,  Range  =  4  km.  Rough  =  0,  Water  Depth  =  100  m 


SNR  (dB) 


Figure  A.20.  BER  vs  SNR  for  16  QAM,  fc  =  8  kHz,  R  =  6.67  kbps,  Range  =  4  km, 

Rough  =  0m,  Water  Depth  =  100m. 


Figure  A.21.  BER  vs  SNR  for  16  QAM,  fc  =  10  kHz,  R  =  8  kbps,  Range  =  2  km, 

Rough  =  0m,  Water  Depth  =  100m. 
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SNR  (dB) 


Figure  A.22.  BER  vs  SNR  for  16  QAM,  fc  =  10  kHz,  R  =  8  kbps,  Range  =  4  km, 

Rough  =  0m,  Water  Depth  =  100m. 


Figure  A. 23.  BER  vs  SNR  for  32  QAM,  fc  =  10  kHz,  R  =  10  kbps,  Range  =  2  km, 

Rough  =  0m,  Water  Depth  =  100m. 
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Figure  A. 24.  BER  vs  SNR  for  32  QAM,  fc  =  10  kHz,  R  =  10  kbps,  Range  =  4  km, 

Rough  =  0m,  Water  Depth  =  100m. 


Figure  A.25.  BER  vs  SNR  for  32  QAM,  fc  =  8  kHz,  R  =  8.33  kbps,  Range  =  2  km, 

Rough  =  0m,  Water  Depth  =  100m. 
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Figure  A.26.  BER  vs  SNR  for  32  QAM,  fc  =  8  kHz,  R  =  8.33  kbps,  Range  =  4  km, 

Rough  =  Om,  Water  Depth  =  100m. 


SNR  =  25dB,  f  =  61cHz,  R  =  5kbps,  Range  =  2  km.  Rough  =  0,  Water  Depth  —  100  m 


SNR  (dB) 


Figure  A. 27.  BER  vs  SNR  for  32  QAM,  fc  =  6  kHz,  R  =  5  kbps,  Range  =  2  km,  Rough 

=  0m,  Water  Depth  =  100m. 
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Figure  A. 28.  BER  vs  SNR  for  32  QAM,  fc  =  6  kHz,  R  =  5  kbps,  Range  =  2  km,  Rough 

=  2m,  Water  Depth  =  100m. 


Figure  A. 29.  BER  vs  SNR  for  32  QAM,  fc  =  6  kHz,  R  =  5  kbps,  Range  =  2  km,  Rough 

=  4m,  Water  Depth  =  100m. 
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Figure  A. 30.  BER  vs  SNR  for  32  QAM,  fc  =  6  kHz,  R  =  5  kbps,  Range  =  4  km,  Rough 

=  0m,  Water  Depth  =  100m. 


Figure  A. 3 1 .  BER  vs  SNR  for  32  QAM,  fc  =  6  kHz,  R  =  5  kbps,  Range  =  4  km,  Rough 

=  2m,  Water  Depth  =  100m. 
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Figure  A. 32.  BER  vs  SNR  for  32  QAM,  fc  =  6  kHz,  R  =  5  kbps,  Range  =  4  km,  Rough 

=  4m,  Water  Depth  =  100m. 


Figure  A. 33.  BER  vs  SNR  for  32  QAM,  fc  =  6  kHz,  R  =  5  kbps,  Range  =  2  km,  Rough 

=  0m,  Water  Depth  =  340m. 
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SNR  =  25dB,  f*c  =  6kHz,  R  =  5kbps,  Range  =  2  km.  Rough  =  2,  Water  Depth  =  340  m 


SNR  (dB) 


Figure  A. 34.  BER  vs  SNR  for  32  QAM,  fc  =  6  kHz,  R  =  5  kbps,  Range  =  2  km,  Rough 

=  2m,  Water  Depth  =  340m. 


Figure  A. 35.  BER  vs  SNR  for  32  QAM,  fc  =  6  kHz,  R  =  5  kbps,  Range  =  4  km,  Rough 

=  0m,  Water  Depth  =  340m. 
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Figure  A. 36.  BER  vs  SNR  for  32  QAM,  fc  =  6  kHz,  R  =  5  kbps,  Range  =  4  km,  Rough 

=  2m,  Water  Depth  =  340m. 


Figure  A. 37.  BER  vs  SNR  for  64  QAM,  fc  =  8  kHz,  R  =  10  kbps,  Range  =  2  km, 

Rough  =  0m,  Water  Depth  =  100m. 


Ill 


Figure  A. 38.  BER  vs  SNR  for  64  QAM,  fc  =  8  kHz,  R  =  10  kbps,  Range  =  4  km, 

Rough  =  0m,  Water  Depth  =  100m. 


Figure  A. 39.  BER  vs  SNR  for  64  QAM,  fc  =  6  kHz,  R  =  5  kbps,  Range  =  2  km,  Rough 

=  0m,  Water  Depth  =  100m. 
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Figure  A. 40.  BER  vs  SNR  for  64  QAM,  fc  =  6  kHz,  R  =  5  kbps,  Range  =  2  km,  Rough 

=  0m,  Water  Depth  =  100m. 


Figure  A. 4 1 .  BER  vs  SNR  for  64  QAM,  fc  =  6  kHz,  R  =  5  kbps,  Range  =  2  km,  Rough 

=  4m,  Water  Depth  =  100m. 
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Figure  A. 42.  BER  vs  SNR  for  64  QAM,  fc  =  6  kHz,  R  =  5  kbps,  Range  =  4  km,  Rough 

=  Om,  Water  Depth  =  100m. 


Figure  A. 43.  BER  vs  SNR  for  64  QAM,  fc  =  6  kHz,  R  =  5  kbps,  Range  =  4  km,  Rough 

=  2m,  Water  Depth  =  100m. 
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Figure  A. 44.  BER  vs  SNR  for  64  QAM,  fc  =  6  kHz,  R  =  5  kbps,  Range  =  4  km,  Rough 

=  4m,  Water  Depth  =  100m. 


Figure  A. 45.  BER  vs  SNR  for  64  QAM,  fc  =  6  kHz,  R  =  5  kbps,  Range  =  2  km,  Rough 

=  0m,  Water  Depth  =  100m. 
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Figure  A. 46.  BER  vs  SNR  for  64  QAM,  fc  =  6  kHz,  R  =  5  kbps,  Range  =  2  km,  Rough 

=  2m,  Water  Depth  =  100m. 


Figure  A. 47.  BER  vs  SNR  for  64  QAM,  fc  =  6  kHz,  R  =  5  kbps,  Range  =  4  km,  Rough 

=  0m,  Water  Depth  =  100m. 


116 


Figure  A. 48.  BER  vs  SNR  for  64  QAM,  fc  =  6  kHz,  R  =  5  kbps,  Range  =  4  km,  Rough 

=  2m,  Water  Depth  =  100m. 
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APPENDIX  B.  MMPE  RESULTS 


Transmission  Loss  (dB  re  lm) 
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Figure  B.  1  Transmission  Loss,  Roughness  0  m,  Range  2  km 

60 

70 

80 

90 

100 

110 

120 

130 

140 

150 


Figure  B.2  Acoustic  Paths,  Roughness  0  m,  Range  2  km ,fc  7  kHz 
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Transmission  Loss  (dB  re  lm) 


6000  6200  6400  6600  6800  7000  7200  7400  7600  7800  8000 


Freq  (Hz) 


Figure  B.3  Transmission  Loss,  Roughness  2  m,  Range  2  km 


Transmission  Loss  (dB  re  lm) 


1.35  Id  1.45  1.5  1.55 

Time  (see) 


Figure  B.4  Acoustic  Paths,  Roughness  2  m,  Range  2  km,  /i  7  kHz 
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Transmission  Loss  (dD  rc  Im) 


6000  6200  6400  6600  6800  7000  7200  7400  7600  7800  8000 


Freq  (Hz) 


Figure  B.5  Transmission  Loss,  Roughness  4  m,  Range  2  km 


Figure  B.6  Acoustic  Paths,  Roughness  4  m,  Range  2  km,  /’.  7  kHz 
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Figure  B.7  Transmission  Loss,  Roughness  0  m,  Range  2  km 


Transmission  Loss  (dB  re  lm) 
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Figure  B.8  Acoustic  Paths,  Roughness  0  m,  Range  2  km,/.  7  kHz 
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Figure  B.9  Transmission  Loss,  Roughness  2  m,  Range  2  km 


Transmission  Loss  (dB  re  lm) 
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Figure  B.  10  Acoustic  Paths,  Roughness  0  m,  Range  2  km ,fc  7  kHz 
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Figure  B.  1 1  Transmission  Loss,  Roughness  0  m,  Range  2  km 


Figure  B.12  Acoustic  Paths,  Roughness  0  m,  Range  2  km,  /'  10  kHz 
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Transmission  Loss  (dB  re  1m) 
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Figure  B.  13  Transmission  Loss,  Roughness  0  m,  Range  2  km 


Figure  B.14  Acoustic  Paths,  Roughness  0  m,  Range  2  km,  /'  12  kHz 
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Figure  B.  15  Transmission  Loss,  Roughness  0  m,  Range  2  km 


Transmission  Loss  (dB  re  Ini) 
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Figure  B.16  Acoustic  Paths,  Roughness  0  m,  Range  2  km,/'  14  kHz 
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APPENDIX  C.  MATLAB  CODE 


A.  TRANSMITTER  CODE 


1.  ofdm  sim  xmitter.m 


%LT  Tiger  Pittman 
%UWA  OFDM  Thesis  m-file 
%Last  modified  4/16/01 
clear  all,  close  all 

fn  =  0;  %running  figure  number  index 

synchdec  =  0; 

%logic  variable  for  synchronization,  0  ==  perfect  synchronization, 
%1  ==  synchronization  using  cyclic  prefix  correlator 


zero  fraction  =  .0250;  %fraction  of  subcarriers  that 

%carry  zeroes,  near  Fs/2  to  avoid  pseudo-aliasing  effects 

[N_bar,  N_subc_act,  N_zerc,  guard_time,  T,  T  fft,  Fs_fft,  zero_fraction,  W,  fc, 
m,  SNR,  delay_spread]  =  parameters(zero_fraction); 

%takes  input  from  user  to  determine  key  parameters  in  the  OFDM  system 

%including  the  number  of  active  subcarriers,  guard  time(correlates  to  cyclic 
prefix  length) 

%OFDM  symbol  period(T,  includes  the  FFT/IFFT  interval  and  the  guard  time), 


noise  =  1 ; 

%logic  variable  for  noise  addition 
noisefactor  =  0.35; 

%if  noise  ==  1  then  AWGN  is  added  with  following  SNR  to  received  signal 
%otherwise  no  noise  is  added 


SNR_awgn  =  SNR*noise_factor; 
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modulate  =  2; 

%if  modulate  =  2  then  double  side  band  modulation  is  performed,  otherwise 
%no  modulation  is  performed 


decision  =  2; 

%detennines  the  infonnation  signal  to  be  used 
%1  is  for  sinusoidal  signals 
%2  is  for  a  converted  voice  wave  file 
%3  is  for  a  ramp  function 


%Signal  generation 

pts  =  2A13;  %number  of  sample  points  (same  as  number  of  quant 

values) 

Fs=250;  %sampling  frequency 

fl  =50;  %tonal  frequency 

f2  =  95; 

al  =  1;  %amplitude 

a2  =  0; 
factor  =  20; 

[info_sig,  fin]  =  signal_generator(pts,Fs,fl,f2, a l,a2, decision,  factor,  fn); 

%calls  to  get  discrete  time  analog  signal 


%Quantization 

bits  =  2A3;  %number  of  bits  used  in  converting  symbol(integer)  to 

binary 

disp([num2str(bits),'  bit  uniform  quantization  is  perfonned  resulting  in  ', 
num2str(2Abits),'  quantization  levels.']) 

[quant  sig,  codebook]  =  quantization(info_sig,  bits); 

%performs  UNIFORM  quantization  of  info  sig  to  create  a 

%discrete  time,  discrete  valued  signal,  quant_sig(i) 
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%Conversion  (quantized  to  binary) 
bitmatrix  =  binary_conversion(quant_sig,  bits); 
%pts  by  bits  matrix  of  binary  data 


%reshape  (binary  matrix  to  vector) 

bit_stream  =  reshape(bit_matrix.',  1,  size(bit_matrix,  1  )*size(bit_matrix,2)); 
%converts  bit  matrix  to  bit  stream  for  DMT  use 

%note:  bit  matrix  is  transposed  since  the  reshape  command  operates  columnwise 
%but  the  data  in  bit  matrix  is  oriented  row  wise 


%QAM  and  Reed  Solomon  parameters 

N  =  2Am  -  1 ;  %Code  word  length  for  Reed_Solomon  symbol  encoding 

q  =  m;  %number  of  bits  per  QAM  symbol,  must  be  <=  m 

Q  =  2Aq;  %QAM  constellation  size,  i.e.  number  of  different 

QAM  symbols 

%Note:  as  long  as  q  is  even,  the  QAM  constellation  is  Gray  coded  (Desired) 

if  (Q  ==  2A6); 

k  =  47;  %RS  message  length  of  Qary  -  QAM  symbols 

elseif  (Q  ==  2A5) 
k  =  23; 

elseif  (Q  ==  2A4) 
k  =  7; 

elseif  (Q  ==  2A3) 
k  =  3; 

elseif  (Q  ==  2A2) 
k=  1; 
else 

error(’QAM  constellation  size  improperly  determined') 
end 

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 

%%%%%%%%%%%%%%%%%%%%%%%%%%%% 

%QAM  Mapping  (binary  word  to  QAM  index) 
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qamcodeddec  =  qam_buffer(bit_stream,  q); 

%buffers  bit  stream  to  q  column  matrix  and  converts 
%to  decimal  format  vector  with  values  0  to  Q-l  of  length 
%ceil(length(bit_stream)/q) 


%Reed  Solomon  encoding 
if  (Q  >=  2A3) 

[rs_code,  zero_symbols]  =  reedsolomon(qam_coded_dec,  N,  k); 
%encodes  QAM  symbols  using  Reed  Solomon  coding 
%rs_code  is  a  matrix  with  each  row  of  length  N 
else 

rscode  =  qamcodeddec; 
zerosymbols  =  0; 
end 


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 


%%%%%%%%%%%%%%%%%%%%%%%%%%%% 


%Data  reshape  (matrix  to  vector) 

code_stream  =  reshape(rs_code,  sizc(rs_codc,  1  )*size(rs_codc,2),  1); 
%reshapes  rs  code  into  column  vector  for  QAM  modulation 
%NOTE:  the  reshape  command  takes  data  columnwise.  However  rs  code 
%is  not  transposed  in  the  reshape  command.  This  effectively 
%interleaves  the  symbols  since  the  subsymbols  are  not  taken 
%in  sequence 


%QAM  Modulation 

[X,  data_carriers,  zero_fdl]  =  qam_modulator(code_stream,  N_bar,  N_subc_act, 
N  zerc,  Q); 

%X  is  the  N(=2*N_bar)  row  by  data_blocks  column  matrix  representing 
%the  spectrum  of  the  transmitted  signal 
Fs_prime  =  2*Fs_fft; 

%sampling  rate  doubles  due  to  the  doubling  of  the  size  of  the  data  which 
%takes  place  when  X  is  created  conjugate  symmetric 
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%IFFT  (converts  complex  QAM  symbols  into  real  time  domain  sequence) 
x  =  inverse_fft(X); 

%N  by  data_blocks  matrix  of  real  time  domain  samples 
%parallel  output  of  the  IFFT 


%Cyclic  extension  application 
L  =  round(guard_time*Fs_primc); 

%L  =  2Aceil(log2(guard_time*Fs_prime));  %expected  length  of  the 

impulse  response  of  the  channel 

%computed  using  the  guard  time  and  the  sampling  rate  at  the  IFFT 

%this  is  the  number  of  smaples  in  the  guard  time,  or  cyclic  extension  period, 
which  must  also  be  an  integer 

x_cyclic  =  cyc_extension(x,  N_bar,  L); 

%inserts  a  cyclic  extension  on  each  of  the  blocks(symbols) 


%Data  reshape  (2*N_bar  by  data_blocks  matrix  to  vector) 
xmit_sig  =  parlel_2_serial(x_cyclic); 

%serial  time  domain  samples,  length  N*data  blocks 


%Modulation  (baseband  signal  modulation  to  carrier  frequency) 
if  (modulate  ==  1) 

[modxmitsig,  Fsxmit,  f_m,  fin]  =  usb_modulator(xmit_sig,  fc,  Fs_prime,  W, 
zerofraction,  fin); 

%modulates  baseband  signal  xmit  sig  to  modulation  frequency,  fc 

%xmit_sig  is  interpolated  to  increase  Fs  to  Fsxmitin  order  to  modulate  w/out 
aliasing 

elseif  (modulate  ==  2) 

[mod  xmit  sig,  Fs  xmit,  f_m,  n_lpf2,  len  xint,  fn]  =  dsb_modulator(xmit_sig,  fc, 
Fs_prime,  W,  zero  fraction,  fn); 

%modulates  the  baseband  signal  using  double  side  band  modulation 
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end 


2.  parameters.m 

function  [N_bar,  N_subc_act,  N_zerc,  guard_time,  T,  T_fft,  Fs_fft,  zero_fraction, 
W,  fc,  m,  SNR,  delay_spread]  =  parameters(zero_fraction) 

%This  function  uses  input  data  to  determine  key  parameters  in  the  OFDM  signal 
generation. 

%Last  modified  7/3/01. 

R  =  1 0 0 0 * ( 5 . 0+ 0/ 3 ) ; % i n p ut( ' W h at  is  the  desired  data  rate  (kbps)?  '); 

W_max  =  1000*7; 
fc  =  1000*6; 
delay_spread  =  le-3*25; 

Pe  =  .0001; 

%one  dimensional  symbol  error  probability,  2  D  symbol  error  probability  =  2*Pe 
SNR_gap3  =  20*logl0(sqrt(2)*erfinv(l-Pe)) 

SNR_gap  =  SNR_gap3  -  10*logl0(3) 

SNR  =  21; 

b_Pe  =  log2(l  +  10A(SNR/10)/10A(SNR_gap/10)); 

%bits  per  subsymbol  that  can  be  carried  at  given  Pe,  SNR_gap,  coding  gain  = 
margin  dB  and  SNR 

guard_time  =  4*delay_spread; 

%guard  time  required  to  prevent  ISI  (sec) 

T  =  6*guard_time; 

%OFDM  symbol  period  including  the  cyclic  prefix 
Rs  =  1/T; 

%OFDM  symbol  rate  (symbols/sec) 
bofdm  =  R*T; 

%bits  per  OFDM  symbol  required  to  get  desired  rate 
N_subc_act  =  ceil(round(b_ofdm)/floor(b_Pe)); 

%number  of  subcarriers  required  for  desired  data  rate  at  achievable  b_Pe 
%these  are  data  carrying  subcarriers 
power2  =  2Aceil(log2(N_subc_act)); 

%determines  the  next  power  of  2  above  N_subc_act 
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if  (((power2-N_subc_act)/power2)  >=  zero_fraction) 

N_bar  =  power2; 
else 

N_bar  =  2Aceil(log2(N_subc_act)  +1); 
end 

N_zerc  =  N_bar  -  N_subc_act; 

%number  of  zero  subcarriers 
zero_fraction  =  N_zerc/N_bar; 

%fraction  of  subcarriers  that  carry  zeroes 
T  fft  =  T  -  guard  time; 

%FFT/IFFT  period 
Fs_fft  =  N_bar/(T_fft); 

%first  guess  at  FFT  sampling  rate  in  FFT  interval 

%based  on  guard  time  =  4*  expected  delay  time  and  an  OFDM  symbol  period  of 
%6*  guard  time 

%this  is  half  the  sampling  rate  required  at  the  DAC  and  LPF  since  the 
%IFFT  doubles  the  number  of  samples  thereby  effectively  interpolating  by  2 
%and  doubling  the  sampling  rate 
nsamp  =  Fs_fft*T; 

%number  of  samples  in  the  OFDM  symbol  period  based  on  Fs  fft 
%must  also  be  an  integer  to  ensure  that  the  subcarriers  are  orthogonal 
if  (isposint(n_samp)~=l) 
while  (isposint(n_samp)~=l) 

Fs_fft  =  Fs_fft  +  2; 
nsamp  =  Fs_fft*T; 
end 
end 

%revises  the  sampling  rate  slightly  in  the  FFT  interval  such  that  there  are  an 
integer 

%number  of  samples  in  the  OFDM  symbol  period 
T_fft  =  N_bar/Fs_fft; 

%revises  the  FFT/IFFT  interval  length  based  on  the  sampling  rate 
%such  that  there  will  be  an  integer  number  of  samples  in  the  interval 
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guard  time  =  T  -  Tfft; 

%revises  the  guard  time  based  on  the  new  T  fft  and  the  known  T 
format  long  e 

SCdelf  =  (T-guard_time)A-l 
%subcarrier  spacing  (Hz) 
format  short  e 

W_mmpe  =  SC_del_f*N_bar 
W  =  SC_del_PN_subc_act; 

%bandwidth  in  Hz 
if  (W  >  W_max) 

disp('Unable  to  meet  requirements  with  given  bandwidth  limitation') 
end 

%system  parameters  output 
m  =  floor(bPe); 
if  (m==5|m==3) 

disp(’QAM  constellation  is  not  Gray  coded’) 
end 

disp('! !!!!!!!!!!  RESULTS  !!!!!!!!!!!!!!!’) 

disp(  ['System  bandwidth  for  DSBSC  transmission  is  ’,num2str(2*W/1000),' 
(kHz)']) 

fimax  =  fc  +  W; 

fimin  =  fc-W; 

if  (fmax>7e3|fmin<200) 

disp('Frequency  spectrum  of  system  is  outside  the  bandwidth  of  the  known 
channel  transfer  function’) 

end 

disp(['System  Transmission  Band  for  USBSC  is  ',  num2str(fmin/1000),'  (kHz)  to  ', 
num2str(fmax/1000),'  (kHz).’]) 

disp(  ['Maximum  Delay  Spread  is  ',  num2str(delay_spread*1000),'  (msec)’]) 

disp(  ['System  uses  ',  num2str(N_subc_act),'  active  subcarriers  with 
’,num2str(N_bar),'  total  subcarriers.']) 

disp(['The  subcarrier  bandwidth  is  ',  num2str(SC_del_f), '  (Hz).']) 
disp(['OFDM  symbol  duration  is  ',  num2str(T*1000),'  (msec).’]) 
disp(['OFDM  Symbol  rate  is  ',  num2str(Rs),'  (symbols/sec).']) 
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disp(  ['Number  of  bits  per  subcarrier  is  ',num2str(m),’.']) 

disp(  [’System  bit  rate  is  num2str(R/1000),  ’  (kbps).']) 

disp(  [’System  Bandwidth  Efficiency  is  ’,num2str(R/(W*2)),'.’]) 

disp(  ['System  carries  ',  num2str(floor(b_Pe)*N_subc_act),’  bits  per  OFDM 
symbol.’]) 


3.  signalgenerator.m 

function  [x,  fin]  =  signal_generator(pts,Fs,fl,f2,al,a2, decision,  factor,  fin) 
%This  function  generates  a  discrete  time  analog  signal. 

%The  signal  is  a  vector  of  length  pts,  from  the  superposition  of 
%2  unit  amplitude  sinusoidal  signals  of  frequency  fl  and  fi2,  sampled 
%at  Fs  Hz. 

%Last  modified  5/15/01 
%global  fin 
if  decision  ==  1 
n=0:pts-l; 

x  =  al *sin(2*pi*fl *n/Fs)+  a2*cos(2*pi*f2*n/Fs);  %information  signal 

spec_info_sig  =  fft(x,2048); 

fn  =  fin+1;  figure(fn),  plot(abs(spec_info_sig)) 

title('FFT  of  sampled  infonnation  signal  (info  sig)') 

elseif  decision  =2 

s  =  wavread(’numnuts'); 

1  =  round(length(s)/factor); 

x  =  s(l:l); 

elseif  decision  ==  3 

x  =  [zeros(  1,100),  0:.2:1,  ones(l,100),  l:-.2:0,  zeros)  1,100)]; 
spec_info_sig  =  fft(x,2048); 

fn  =  fin+1;  figure(fn),  plot(abs(spec_info_sig)) 
title(’FFT  of  sampled  infonnation  signal  (info  sig)') 
end 
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4.  quantization. in 

function  [symbol,  codebook]  =  quantization/infosig,  bits) 

%Last  modified  3/19/01 

M  =  2Abits;  %number  of  uniform  quantization  levels 
%must  determine  the  method  for  partitioning  to  use  in  general 
maxsample  =  max(abs(info_sig)); 

%absolute  maximum  value  of  input  signal 

part=-/M-l)/M*max_sample:2/M*max_sample:(M-l)/M*.9*max_sample; 

%partitions 

codebook  =  -((M-l)*2  +  l)/(2*M)*max_sample:2/M*max_sample:((M- 

l)*2+l)/(2*M)*max_sample; 

symbol  =  quantiz/infosig,  part); 

5.  dec_2_bin.m 

function  bit  stream  =  dec_2_bin(qam_sym,  bits,  q) 

%Last  modified  3/19/01 
bit_matrix  =  de2bi(qam_sym',q); 

%binary  matrix  length(qam_sym)  by  q 

bit_stream_long  =  reshape(bit_matrix,  length(qam_sym)*q,  1); 
num_zero_bits  =  length/ qam_sym)*q  -  floor((length(qam_sym)*q)/bits)*bits; 
%number  of  zero  bits  inserted  in  xmitter 

bitstream  =  bit_stream_long(num_zero_bits+ 1  :length(bit_stream_long)); 

6.  qam  buffer.m 

function  qamdec  =  qam_buffer(bit_stream,  q) 

%This  function  buffers  the  bit  stream  into  a  matrix  form  for  QAM  encoding. 
%Last  modified  3/19/01 

add_length  =  q*(ceil(length(bit_stream)/q))-length(bit_stream); 

%detennines  the  number  of  zeroes  needed  to  be  added  to  bit  stream 
%such  that  buffered  matrix  will  be  rectangular  and  include  all  of  the  data. 
bit_stream_long  =  [zeros/ 1  ,add_lcngth),bit_strcam]; 

%bit_stream  with  zeroes  added  to  left  end  i.e.  at  1,2,3... 
buffer_mat  =  reshape(bit_stream_long,length(bit_stream_long)/q,q); 
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%buffers  bit  stream  into  q  column  matrix  for  QAM  modulation 
%data  is  reshaped  columnwise 

qamdec  =  bi2de(buffer_mat);  %converts  buffered  binary  data  to 
corresponding 

%decimal  integer  where  the  integers  are  from  0  to  Q-l(2Aq  -  1) 

%this  reduces  the  sampling  frequency  since  there  are  now  less  samples 
%least  sig  bit  assumed  in  column  1 


7.  reedsolomon.m 

function  [code,  zerosymbols]  =  reedsolomon(qam_coded_dec,  N,  k) 

%This  function  performs  Reed  Solomon  encoding  of  the  decimal  format 
%Gray  coded  Qary  QAM  symbols  in  qamcodeddec 
%k  is  the  message  length  in  symbols,  N  is  the  code  word  length  in  symbols 
%Last  modified  3/19/01 

fullblocks  =  floor(length(qam_coded_dec)/k); 

%number  k  length  blocks  filled  with  data 

zero  symbols  =  (full  blocks  +  l)*k  -  length(qam_coded_dec); 

%number  of  zero  blocks  needed  to  fill  final  block 
qamcodeddeclong  =  [zeros(zero_symbols,l);qam_coded_dec]; 

%QAM  code  with  zeroes  at  top  of  column  vector 
qamkmat  =  reshape(qam_coded_dec_long,full_blocks+l,k); 

%full_blocks  +  1  row  by  k  column  QAM  symbol  matrix  for  block  coding 
code  =  encode(qam_k_mat,  N,  k,  'rs/decimal'); 

%full_blocks  +  1  row  by  N  column  matrix  of  Qary-QAM  encoded  symbols 
%RS  encoded 

8.  qammodulator.m 

function  [X,  data  carriers,  zero  fill]  =  qam_modulator(code_stream,  N  bar, 
N_subc_act,  N_zerc,  Q) 

%Last  modified  5/15/01 

[inphase,  quad]  =  qaskenco(code_stream,Q); 

%encodes  integer  decimals  into  quadriture  and  inphase  components  of  Q-ary 
QAM 
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[phaseqam,  magqam]  =  cart2pol(inphase,quad); 

%converts  QAM  subsymbols  to  magnitude  and  phase(radians)  for  IFFT 
zero_carriers  =  N_zerc;  %ceil(zero_fraction*N_bar); 

%number  of  zero  subcarriers  near  Fs/2 
data_carriers  =  N_bar  -  zero_carriers; 

%number  of  subcarriers  with  encoded  data 
data_blocks  =  ceil(length(phase_qam)/data_carriers); 

%number  of  blocks(  i.e.  symbols)  neeeded  for  given  N  bar  and  zero  fraction 
zerofill  =  data_blocks*data_carriers  -  length(phase  qam); 
colz  =  floor(zero_fill/data_blocks); 

%number  of  zeroes  to  be  added  into  QAM  symbol  matrix  for  columns  2  : 
datablocks 

collz  =  zero  fill  -  (data_blocks-l)  *  col  z; 

%number  of  zeroes  to  be  added  to  the  first  column  in  the  QAM  symbol  matrix 

zero  fill  =  [col  lz,  ones(l,  data_blocks+2-l)*col_z]; 

phase  cl  =  phase_qam(l:(data_carriers  -  col  lz)); 

phase  col  =  phase_qam(l+(data_carriers  -  col_lz):length(phase_qam)); 

col_phase_mat  =  reshape(phase_col,  datacarriers  -  col  z,  data_blocks-l); 

magcl  =  mag_qam(l:(data_carriers  -  col  lz)); 

magcol  =  mag_qam(l+(data_carriers  -  col_lz):length(mag_qam)); 

col  mag  mat  =  reshape(mag_col,  data  carriers  -  col  z,  data_blocks-l); 

tpcl  =  [phase_cl(l:col_lz).';zeros(l,col_lz)]; 

tpc2  =  reshape(tpc  1 ,  2*col_lz,  1); 

phaseclO  =  [tpc2;  phase_cl(col_lz+l:length(phase_cl))]; 

tmcl  =  [mag_cl(l:col_lz).';zeros(l,col_lz)]; 

tmc2  =  reshape(tmc  1 ,  2*col_lz,  1); 

mage  10  =  [tmc2;  mag_cl(col_lz+l:length(mag_cl))]; 

for  col  =  1  :  data  blocks  -  1 

tpcola  =  [col_phase_mat(l:col_z,  col).';  zeros(l,  col  z)]; 
tpcolb  =  reshape(tpcola,  2*col_z,  1); 

phasecolO(:,col)  =  [tpcolb;col_phase_mat(col_z+l:size(col_phase_mat,l),  col)]; 
tmcola=  [col_mag_mat(l:col_z,  col).';  zeros(l,  col  z)]; 
tmcolb  =  reshape(tmcola,  2*col_z,  1); 
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magcolO(:,col)  =  [tmcolb;col_mag_mat(col_z+l:size(col_mag_mat,l),  col)]; 
end 

phaseqammat  =  [phaseclO,  phasecolO]; 
magqammat  =  [mage  10,  magcolO]; 
for  k  =  1  :data_blocks 

X_half(:,k)  =  zeros(N_bar+l,l); 
for  i=l  :data_carriers; 

X_half(i+l,k)  =  mag_qam_mat(i,k).*cxp(].*phasc_qam_mat(i,k)); 
end 

X_half(l,k)  =  real(X_half(data_carriers+l,k)); 
X_half(data_carriers+l,k)  =  imag(X_half(data_carriers+l,k)); 

X(:,k)  =  [X_half(:,k);  conj'(flipud(X_half(2:N_bar,k)))]; 
end 

count  =  0; 

for  r  =  data_carriers+ 1 :  - 1 : 1 
for  c  =  size(X_half,2):-l :  1 
count  =  count  +  1 ; 

X_pilot_half(count,l)  =  X_half(r,c); 
end 
end 

X_half(:,data_blocks+l)  =  zeros(N_bar+l,l); 

X_half(:,data_blocks+2)  =  zeros(N_bar+l,l); 

X_half(2:data_carriers+l,data_blocks+l)  =  X_pilot_half(l  :  data_carriers); 

X_half(  1  ,data_blocks+ 1 )  =  real(X_half(data_carriers+ 1  ,data_blocks+ 1 )); 

X_half(data_carriers+ 1  ,data_blocks+ 1 ) 
imag(X_half(data_carriers+ 1  ,data_blocks+ 1 )); 

X_half(  1  :data_carriers+ 1  ,data_blocks+2) 
flipud(X_half(  1  :data_carriers+ 1  ,data_blocks+l)); 

for  c  =  1 :2 

for  r  =  1  :data_carriers+l 
if  X_half(r,2)  ==  0 
X_half(r,data_blocks+c)  =  0; 
end 
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end 

end 

for  k  =  data_blocks+l:data_blocks+2 
X(:,k)  =  [X_half(:,k);  conj'(flipud(X_half(2:N_bar,k)))]; 
end 

X_pilots  =  X(:,  data_blocks+l:size(X,2)); 

X  =  [X_pilots,  X(:,l:data_blocks)]; 
for  k  =  1  :data_blocks+2 
for  i=l  :data_carriers; 
if  (k  <=  datablocks) 

X_half(i+l,k)  =  mag_qam_mat(i,k).*exp(j.*phase_qam_mat(i,k)); 

%creates  first  half  of  complex  symbol,  other  half  is  the 

conjugate 

%symmetric  counterpart 
else 

X_half(i+l,k)  =  1; 
end 
end 

X_half(l,k)  =  real(X_half(data_carriers+l,k)); 

%first  channel  is  the  real  part  of  the  last  symbol 
if  k  <=  datablocks 

X_half(data_carriers+l,k)  =  imag(X_half(data_carriers+l,k)); 
end 

X(:,k)  =  [X_half(:,k);  conj(flipud(X_half(2:N_bar,k)))]; 
end 


9.  inverse  fft.m 

function  x_n  =  inverse  fft(X) 

%Last  modified  4/13/01 
x_n  =  ifft(X); 

%N  time  domain  samples,  real  valued  from  symmetric  fireq  spectrum,  X 
%x_n  is  a  matrix  the  same  size  as  X,  representing  the  time  domain  samples 
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%to  be  transmitted  for  each  block 
max_imag  =  max(max(imag(x_n))); 
if  (max_imag>=le-14) 

error(’Time  domain  output  from  IFFT  not  real  to  machine  accuracy') 
end 


10.  cycextension.m 

function  x_cycex  =  cyc_extension(x,  N_bar,L) 
%Last  modified  3/19/01 
for  i  =  1  :size(x,2) 

x_cycex(:,i)  =  [x(2*N_bar-L+l:2*N_bar,i);x(:,i)]; 
end 


11.  parlel_2_serial.m 

function  x_ser  =  parlel_2_serial(x_cyclic) 

%This  function  converts  the  parallel  output  of  the  IFFT  to  serial  form. 
%x_ser(l)  represents  the  last  sample  to  be  transmitted 
%Last  modified  3/19/01 

x  =  reshape(x_cyclic,  1,  size(x_cyclic,l)*size(x_cyclic,2)); 
x_ser  =  real(x); 

%real  command  necessary  to  correct  for  small  imaginary 
%components  in  x  from  machine  error 

%converts  symbol  from  parallel  to  serial  w/  x(N,data_blocks)leading  the  signal 
%and  x(l,l)  trailing 


12.  dsbmodulator.m 

function  [mod  sig,  Fsmod,  f_m,  n_lpf2,  lenxint,  fn]  =  dsb_modulator(xmit_sig, 
fc,  Fs_prime,  W,  zerofraction,  fn); 

%This  function  modulates  the  baseband  signal,  xmitsig  to  the  carrier  frequency, 
fc,  using 

%double  sideband  suppressed  carrier  (DSBSC)  modulation. 
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%The  baseband  signal  is  interpolated  as  necessary  to  meet  Nyquist  frequency 
requirements 

%Last  modified  4/5/01 


%plots 

N=8 1 92;  %fft  length 

f_p=(0:N-l)*Fs_prime/N;  %index  for  plotting  fft's 


%ERROR  CHECK  TO  ENSURE  NO  ALIASING 


if  (fc  <  W) 

error(Aliasing  will  occur  due  to  improper  fc  and  W  combination,  fc  must  be 
greater  than  W.') 

end 


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 


%%%%%%%%%%%%% 


%  UPSAMPLING  TO  INCREASE  THE  SAMPLING  RATE 

K  =  1.0; 

%  margin  above  the  Nyquist  rate 
I  =  ceil(2*K*(fc+W)/Fs_prime); 

/^interpolation  constant 

[xint,  fin]  =  interp_mod(xmit_sig,  I,  fin); 

%interpolated  xmit  signal,  baseband,  DSB 
lenxint  =  length(xint); 

Fsmod  =  I*Fs_prime; 

%  sampling  frequency  of  xmitted  and  modulated  signal 
f_m  =  (0:N- 1  )*Fs_mod/N; 

%frequency  index  for  interpolated  signal  spectrum 


%POST  INTERPOLATION  LPF 
fracpass2  =1.2; 

%  fraction  of  W  for  pass  band 
fracstop2  =  1.375; 

%  fraction  of  W  for  stop  band 
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digf_pass2  =  2*pi*fracpass2*W/(Fs_mod);  %  relative  pass  frequency 
digf_stop2  =  2*pi*fracstop2*W/(Fs_mod);  %  relative  stop  frequency 
xtion2  =  digf_stop2-digf_pass2;  %transition  region  width 

fpass2  =  fracpass2*W/(.5*Fs_mod); 

%fraction  of  Fs/2  for  cutoff  frequency  in  FIR  filter  from  firl.m 

n_lp!2  =  round(  8*pi/xtion2  );  %hamming  filter  order 

if  (isposint(n_lpl2/2)  ~=  1 )  %ensures  filter  has  odd  number  of 

coefficients 

n_lpf2  =  n_lp(2+ 1 ; 

end 

[Nh2,Dh2]  =  firl(  n_lpf2  ,  fpass2  );  %Hamming  LP  filter  numerator  and 

denominator  coefficients 

xintlpf  =  conv(Nh2,  x_int); 


mod  =  cos(2*pi*fc/Fs_mod*[l:length(x_int_lpf)]); 
modsig  =  mod.*x_int_lpf; 

mod_sig_spec  =  fft(mod_sig,  N);%,  2Aceil(log2(length(mod_sig)))); 
xmit_sig_spec  =  fft(xmit_sig,  N);%,  2Aceil(log2(length(xmit_sig)))); 


B.  RECEIVER  CODE 


1.  ofdm  sim  receiver.m 

%LT  Tiger  Pittman 
%OFDM  Simulation:  receiver 
%Last  Modified  5/16/01 
depth  =  0; 

%depth  =  0  corresponds  to  shallow  water  model,  100m  depth 
%depth  =  1  corresponds  to  deep  water  model,  340m  depth 
range  =  1 ; 

%range  0  corresponds  to  2  km  transmission  range,  1  to  4  km  range 
rough  =  0; 
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%can  be  0,2  or  4 

%0  for  all  allowed  fc,  2  for  deep  water  and  fc  =  7  only, 
%2  and  4  for  shallow  water  at  fc  =  7  only 
recvd_sig  =  mod_xmit_sig; 


%  ADDITIVE  WHITE  GAUSSIAN  NOISE 
if  noise  ==  1 

%Ps=(sum(recvd_sig(  1 : 1 000).  A2)/ 1 000); 

Ps=(sum(recvd_sig(l:length(recvd_sig)).A2)/length(recvd_sig)); 

Pn  =  Ps/((10A(SNR_awgn/10))); 

awgn  =  randn(l,length(recvd_sig))*sqrt(Pn); 

recvd_sig  =  recvd_sig  +  awgn; 

%  verification  of  SNR  in  time  domain 
Pnoise=(sum(awgn(  1 : 1000).A2)/1000); 
SNR_est=10*logl0(Ps/Pnoise);  %  checks  SNR 
display(SNR_est) 
end 


%DEMODULATION,  (includes  low  pass  filtering  and  decimation) 
if  (modulate  ==  1) 

[recvd_sig,  fin]  =  rcvr_usb_demod(recvd_sig,  Fs_xmit,  fc,  mod_xmit_sig, 
W,  Fs_prime,  fin); 

elseif  (modulate  ==  2) 

[recvd_sig,  fin]  =  rcvr_dsb_demod(recvd_sig,  Fs_xmit,  fc,  mod_xmit_sig,  W, 
Fs_prime,  n_lpf2,  len_xint,  synchdec,  fin); 

end 


%performs  time  and  frequency  synchronization 
blocklength  =  size(x_cyclic,l); 
num_blocks  =  size(X,2); 
if  synchdec  ==  1 
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[sync_sig,  synch,  fn]  =  synchronizer(recvd_sig,  block_length,  num_blocks, 
N_bar,  L,  SNR_awgn,  fn); 

else 

sync_sig  =  recvd_sig; 
end 


%Serial  to  parallel  conversion 

par_sig_r  =  ser_2_parlel(sync_sig,  L,  N_bar); 

%parallel  blocks(symbols)  of  received  data  with  cyclic  extension 
%still  attached  at  top  of  each  symbol(column) 


%Cyclic  extension  removal 

par  =  par_sig_r(L  +  1 :2*N_bar  +  L,  :); 


%FFT  of  received  signal 
X_r_eb  =  fft(par); 

%FFT  of  received  symbols,  includes  channel  equalization  blocks 


%Channel  application 

[X_r_ebch,  rcvr_depth,  fn]  =  mmpe_channel(X_r_eb,  N_bar,  N_subc_act,  fc, 
depth,  range,  rough,  fn); 


%performs  channel  equalization 

X_r  =  equalization(X_r_ebch,  X,  N_subc_act,  N_bar,  fn); 


%QAM  Demodulator 

[code,  inphase  r,  quad  r,  fn]  =  qam_demodulator(X_r,  N  bar,  zero  fraction,  Q, 
data_carriers,  N,  zero  fill,  fn); 

%demodulates  X_r  to  give  Reed  Solomon  encoded  integers,  0:Q-1 
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code_deint  =  deinterleaver(code,  N); 
%deinterleaves  the  received  symbols 


if  (Q  >=  2A3) 

qam_sym  =  rsdecoder(code_deint,  N,  k,  zero_symbols); 
%decodes  RS  encoded  symbols 

else 

qamsym  =  codedeint; 
end 


bit_stream_r  =  dec_2_bin(qam_sym,  bits,  q); 

%reshapes  qam  sym  to  q  column  matrix  for  conversion  to  binary 
%removes  zero  bits  inserted  in  xmitter 


buffer_r  =  reshape(bit_stream_r,  bits,  length(bit_stream_r)/bits)'; 
%reshapes  bitstreamr  to  bits  column  matrix  for  conversion  to 
%decimal  integers  for  signal  reconstruction 


quantint  =  bi2de(buffer_r); 

%converts  binary  words  of  length  bits  to  decimal  integer 


infosigr  =  codebook(quant_int  +  1); 


signals2(info_sig,  info_sig_r,  code,  code_stream,  Q,  X_r,  fin) 
BER  =  error_check(bit_stream,  bit  stream  r); 
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2.  rcvr  dsb  demodulator.m 


function  [decfrsig,  fn]  =  rcvr_dsb_demod(recvd_sig,  Fsxmit,  fc, 
mod_xmit_sig,  W,  Fs_prime,  n_lpf2,  len_xint,  synchdec,  fn); 

%This  function  demodulates  the  DSBSC  signal  from  the  channel  transfer  function 

%Last  modified  5/9/01 


%DEMODULATION 

demodsig  =  cos(2*pi*fc*(l:length(recvd_sig))/Fs_xmit); 
%demodulating  signal  to  combine  upper  side  band  spectrums  to  recreate 
%baseband  signal 

demod_rec_sig  =  recvd_sig.*demodsig; 

%demodulated  time  domain  signal 


N  =  2*8192;  %length(demod_rec_sig); 
f  =  ( 0 : N - 1  )*Fs_xmit/N;  %index  for  plotting  fft's 

%frequency  index  of  received  signal  spectrum 


%LOW  PASS  FILTERING 


fracpass  =  1.075; 

%  fraction  of  W  for  pass  band 
fracstop  =  1.2; 

%  fraction  of  W  for  stop  band 

digf_pass=2  *pi*  fracpass*  W/(Fs_xmit);  %  relative  pass  frequency 

dig  f_stop=2*pi*  fracstop*  W/(Fs_xm  it);  %  relative  stop  frequency 

if  (digf_stop>=.99*pi) 

error('Stop  band  exceeds  .99  of  the  Sampling  frequency') 
end 

xtion=digf_stop-digf_pass;  %transition  region  width 

fpass  =  fracpass*  W/(. 5  *Fs_xmit); 

%fraction  of  Fs/2  for  cutoff  frequency  in  FIR  filter  from  firl.m 
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n_lpf  =  round(  8*pi/xtion  );  %hamming  filter  order 

if  (isposint(n_lpf/2)  ~=  1)  %ensures  filter  has  odd  number  of 

coefficients 

nlpf  =  n_lpf+l; 

end 

[Nh,Dh]  =  firl(  n  lpf,  fpass);  %Hamming  LP  filter  numerator  and 

denominator  coefficients 

filt  rsig  =  conv(Nh,  demod_rec_sig); 

if  synchdec  ==  0 

timedelaymod  =  n_lpf2; 

timedelaydemod  =  nlpf; 

time_delay_tot  =  floor(time_delay_mod/2)  +  floor(time_delay_demod/2); 

filt_rsig  =  filt_rsig(l,  time_delay_tot+l:len_xint+time_delay_tot); 

%THIS  COMMAND  PERFORMS  PERFECT  SYNCHRONIZATION 

end 


%DECIMATION 


D  =  round(Fs_xmit/Fs_prime); 
if  (isposint(D)~=l) 

error('Decimation  factor  not  an  integer') 
end 

if  (isposint(length(  f11trsig)/D)~=  1 ) 
if  D  ==  2 

filt_rsig  =  [filt  rsig,  0];  %appends  zero  to  signal  so  that  it  will  decimate 
evenly 

elseif  D  >=  3 

lenl  =  length(filt_rsig)  +  1; 
len2  =  length(filt_rsig)  +  2; 
len3  =  length(filt_rsig)  +  3; 
if  isposint(lenl/D) 
filt  rsig  =  [filt_rsig,  0]; 
elseif  isposint(len2/D) 
filt  rsig  =  [filt  rsig,  0  0]; 
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else 

filt  rsig  =  [filt  rsig,  0  0  0]; 
end 
end 
end 

if  D>=5 

error(’Interpolation/Decimation  factors  greater  than  4  requires  code  adjustment 
in  rcvrdsbdemod.m') 

end 

dec_mat  =  reshape(fdt_rsig,  D,  length(fdt_rsig)/D); 
dec_f_rsig  =  dec_mat(l,:); 

%decimated  received  signal,  returns  the  sampling  frequency  to  Fs_prime  = 
Fs_xmit/D 

f_p  =  (0:N-l)*Fs_prime/N;  %index  for  plotting  fft's 

%frequency  index  of  received  signal  spectrum 


3.  synchronizer.m 


function  [sync_sig,  synch,  fn]  =  synchronizer(recvd_sig,  block_length, 
num_blocks,  N_bar,  L,  SNR_awgn,  fn) 

%Correlator  synchronizer 

%last  modified  6/21/2001 


%  ADDITIVE  WHITE  GAUSSIAN  NOISE 


%Ps=(sum(recvd_sig(  1 : 1 000) .  A2)/ 1 000); 
Ps=(sum(recvd_sig(l:length(recvd_sig)).A2)/length(recvd_sig)); 
%n=0: 1  :N-1; 

%S=sin(2*pi*  16.6*n/FS)+sin(2*pi*  17.9*n/FS); 

Pn  =  Ps/((10A(SNR_awgn/10))); 

awgn  =  randn(l,length(recvd_sig))*sqrt(Pn); 

%  verification  of  SNR  in  time  domain 
Pnoise=(sum(awgn(  1 : 1 000).  A2)/ 1 000); 
SNR_est_synch=10*logl0(Ps/Pnoise);  %  checks  SNR 
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display(SNR_est_synch) 
lengnoisel  =  1000; 
leng_noise2  =  1500; 

r  =  [sqrt(Pn)*randn(l,  leng_noisel),  recvd_sig,  sqrt(Pn)*randn(l,  leng_noise2)  ]; 


for  k  =  1  :length(r)-2*N_bar-2*L+ 1  %loops  over  entire  length  of  data 

r_sig  =  r(l,  k:k+L-l); 

r_del_conj_sig  =  conj(r(l,  k+2*N_bar:k+2*N_bar+L-l)); 
corr_sig(k)  =  max(xcorr(r_sig,  r_del_conj_sig)); 
end 

[synch,  fin]  =  peakfmder(corr_sig,  fin); 

%last  block  identifies  the  ofdm  block  length 

%previous  blocks  identify  the  discrete  time  prior  to  the  beginning  of  the  data 
blocks 

%i.e.  synch(l)  is  the  beginning  of  the  first  data  block  and  synch(end)  is  the  ofdm 
block  length 

%performance  check 

block_length_error  =  abs(block_length  -  synch(length(synch))); 
peak  count  error  =  numblocks  -  (length( synch)  -  1); 
block  start  error  =  abs(leng_noisel  +  1  -  synch(l)); 

disp(  ['Block  length  error  =  num2str(block_length_error),’  discrete  time  units.']) 

disp(  ['Peaks  missed  =  ',  num2str(peak_count_error),’.']) 

sync_sig  =  r(synch(l)+l  :  synch(l)  +  1  +  num_blocks*block_length  -  1); 


a.  peakfinder.m 

function  [synch,  fin]  =  peakfinder(corr_sig,  fin) 

%last  modified  5/16/01 
peakindex  =  0; 

%peaks  =  find(corr_sig>mn  +  2.5  *sd); 
fracnum  =  0; 

tl  =  .4;  t2  =  .6;  dt  =  .05;  %threshholds  and  differential 

for  fraction  =  1 1  :dt:t2 
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fracnum  =  fracnum  +  1 ; 

%fraction  =  0.7; 

peaks  =  find(corr_sig  >  fraction*max(corr_sig)); 
count  =  0; 
peaknum  =  0; 

while  count  <  length(peaks) 
samepeak  =  1 ; 

peaknum  =  peaknum  +  1 ;  %index  for  number  of  peaks  found 
lookstart  =  count  +  1 ; 
while  same_peak  ==  1 ; 
count  =  count  +  1 ; 
if  (count  >=  length(peaks)) 
break 
end 

if  (peaks(l+count)  -  peaks(count)  ==  1) 
same_peak  =  1 ; 
else 

samepeak  =  0; 
end 
end 

%while  and  if  loop  determine  the  indices  of  corr  sig 
corresponding  to  the  1st  peak  which  is  where 

%the  first  OFDM  symbol  begins 

%count-l  after  while  loop  is  the  last  index  of  the  1st  peak, 
nextpeakstart  =  count  +  1 ; 

[mag_peak(peaknum),  peak(peaknum)]  = 

max(corr_sig(peaks(lookstart):peaks(count))); 

peak_index(fracnum,  peaknum)  =  peak(peaknum)  + 

peaks(lookstart)  -  1 ; 

end 

end 

%disp(fraction) 

%disp(peak_index) 
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numtruepeaks  =  0; 
for  1  =  1  :size(peak_index,2) 
if  peak_index(l,l)  ~=  0 

matches  =  fmd(peak_index(l,l)  == 

peak_index(2:size(peak_index,l),:)); 

if  length(matches)  ==  size(peak_index,l)  -  1 

numtruepeaks  =  numtruepeaks  +  1 ; 

truepeakindex(numtruepeaks)  =  peak_index(l,l); 

end 

end 

end 

%break 

blocklength  =  truepeakindex(2)  -  truepeakindex(l)  -  1; 
synch  =  [truepeakindex,  block  length]; 

fn=fn+l;  figure(fn),  plot(corr  sig),  title(’Correlator  Synchronizer  Output'), 
hold  on 

plot(truepeakindex,  corrsig(truepeakindex),  'r*') 
fork  =  tl:dt:t2 

line([l,  length(corr_sig)],  [k*max(corr_sig),  k*max(corr_sig)]) 
end 

hold  off 

4.  mmpe  channel.m 

function  [X_ch,  rcvr_depth,  fin]  =  mmpe_channel(X_r_eb,  N_bar,  N_subc_act,  fc, 
depth,  range,  rough,  fn) 

%This  function  applies  the  transfer  function  from  the  MMPE  model  to  the 
received 

%QAM  symbols.  X_r_ebch(:,i)  =  X_r_eb(:,i).*H(z),  where  i  is  the  OFDM  block 
index 

% 

%Last  modified  6/4/01 
if  N_bar  ==  1024 
if  fc  ~=  6e3 

error('Carrier  frequency  does  not  correspond  to  model  carrier  frequency') 
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end 

if  (depth  ==  0  &  range  ==  0  &  rough  ==  0) 
load  f7r02k.mat; 

elseif  (depth  ==  0  &  range  ==  0  &  rough  ==  2) 
load  f7r22k.mat; 

elseif  (depth  ==  0  &  range  ==  0  &  rough  ==  4) 
load  f7r42k.mat; 

elseif  (depth  ==  0  &  range  ==  1  &  rough  ==  0) 
load  f7r04k.mat; 

elseif  (depth  ==  0  &  range  ==  1  &  rough  ==  2) 
load  f7r24k.mat; 

elseif  (depth  ==  0  &  range  ==  1  &  rough  ==  4) 
load  f7r44k.mat; 

elseif  (depth  ==  1  &  range  ==  0  &  rough  ==  0) 
load  f7r02kd.mat; 

elseif  (depth  ==  1  &  range  ==  0  &  rough  ==  2) 
load  f7r22kd.mat; 

elseif  (depth  ==  1  &  range  ==  1  &  rough  ==  0) 
load  f7r04kd.mat; 

elseif  (depth  ==  1  &  range  ==  1  &  rough  ==  2) 
load  f7r24kd.mat; 

end 

elseif  N  bar  ==  2048 


if  (fc  ==  8e3  &  range  ==  0) 
load  fl0r02k.mat; 
elseif  (fc  ==  8e3  &  range  ==  1) 
load  fl0r04k.mat; 
elseif  (fc  ==  10e3  &  range  ==  0) 
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load  fl2r02k.mat; 
elseif  (fc  ==  10e3  &  range  ==  1) 
load  fl2r04k.mat; 
elseif  (fc  ==  12e3  &  range  ==  0) 
load  fl4r02k.mat; 
elseif  (fc  ==  12e3  &  range  ==  1) 
load  fl4r04k.mat; 
end 
end 

r  =  size(press,l);  %number  of  depths  in  model,  first  half  are  in 

acoustic  channel,  i.e.  water 

[v  dep_index]  =  max(sum(abs(press(l:r/2,:)),2));  %dep_index  indetifies 

the  depth  with  the  greatest 

%return  over  all  frequencies 

H  =  press(dep_index,:);  %places  the  receiver  at  the  optimal  depth 

%H  =  press(10,:); 

Hhalf  =  zeros(length(H)+l,l); 
fork=  l:length(H) 

H_half(k+1)  =  H(k); 
end 

H_half(l)  =  real(H_half(length(H_half))); 

H_half(length(H_half))  =  imag(H_half(length(H_half))); 

H  full  =  [H  half;  conj(flipud(H_half(2:N_bar)))]; 

%transfer  function  to  be  applied  to  received  QAM  symbols 
h  =  flipud(ifft(H_full)); 

rcvr_depth  =  ( 1  -depth)*(dep_index/(r/2))*  1 00  +  depth*(dep_index/(r/2))*340; 
%receiver  depth  in  meters 
for  k  =  1  :size(X_r_eb,2) 

X_ch(:,k)  =  X_r_eb(:,k).*H_full; 

End 

5.  equalization.m 

function  X_r  =  equalization(X_r_eb,  X,  N_subc_act,  N_bar,  fn) 
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%performs  channel  equalization  using  pilot  blocks 
%last  modified  5/15/01 

H_est_mat_real  =  X_r_eb(l:N_subc_act+l,  1:2)./  ... 

X(  1  :N_subc_act+ 1 ,  1:2); 

H_est_mat_imag  =  X_r_eb(size(X_r_eb,l)  -  N_subc_act  +  1:  size(X_r_eb,l), 

1:2)./... 

X(size(X,l)  -  N_subc_act  +  1:  size(X,l),  1:2); 

%estimate  of  channel  freqency  response  from  2  unit  magnitude  OFDM  blocks 

%estimate  is  the  average  of  the  2  columns  and  is  used  for  channel  equalization 

%real  is  estimate  for  first  N_subc_act  subcarriers,  and  imag  is  for  last  N_subc_act 
subcarriers 

%equalization  is  only  performed  over  the  active  subcarriers 
Hestreal  =  mean(H_est_mat_real.').'; 

Hestimag  =  mean(H_est_mat_imag.').'; 

%2*N_subc_act  column  vector  representing  the  channel  estimate  over  the  actice 
subcarriers 

X_r  =  X_r_eb(:,  3:size(X_r_eb,2)); 

%removes  channel  estimation  blocks  assuming  the  first  2  blocks  are  the  pilot 
blocks  and  no  others 

forkl  =  l:size(X_r,2); 

X_r_real(:,kl)  =  X_r(l:N_subc_act+l,kl)./H_est_real; 

X_r_imag(:,kl)  =  X_r(size(X_r,l)  -  N_subc_act  +  1: 

size(X_r,  l),kl)./H_est_imag; 

X_r(:,kl)  =  [X_r_real(:,kl);  zeros(2*N_bar  -  2*N_subc_act  -  1,1); 

X_r_imag(:,kl)]; 

end 

%equalizes  the  data  for  the  channel 


6.  qamdemodulator.m 

function  [code,  inphaser,  quadr,  fn]  =  qam_demodulator(X_r,  Nbar, 
zero  fraction,  Q,  data  carriers,  N,  zero  fill,  fn) 

%Last  modified  5/11/01 

X_half_r  =  X_r(l:data_carriers+1,:); 

%strips  conjugate  symmetric  and  zero  subcarriers  from  data 
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for  k  =  1 :  size(X_half_r,2) 
for  i  =  1  :data_carriers 
mag_r(i,k)  =  abs  (X_half_r(i+l,k)); 

%determines  the  magnitude  of  the  QAM  symbol 
phase_r(i,k)  =  angle(X_half_r(i+l,k)); 

%determines  the  phase  of  the  QAM  symbol 
end 

mag_r(data_carriers,k)  =  sqrt(X_half_r(l,k)A2  + 

real(X_half_r(data_carriers+ 1  ,k))A2); 

%combines  the  Re(N_bar+l)  and  the  Imag(N_bar+l)  that  were  separated 

%in  the  qam  modulator  to  get  the  magnitude  if  the  final  symbol 

phase_r(data_carriers,k)  =  atan2(real(X_half_r(data_carriers+ 1  ,k)), 

X_half_r(  1  ,k)); 

%combines  the  Re  and  Im  at  N  bar  +  1  to  get  the  phase  for  the  final 

symbol 

%real  command  necessary  to  account  for  small  imaginary  component 
resulting  from 

%machine  error 

end 

%result  of  these  2  loops  is  a  data  carrier  row  matrix  representing  the  QAM 
modulated 

%symbols  where  each  column  is  the  kth  symbol  and  each  row  is  the  ith 
subsymbol 


tmagl(:,l)  =  mag_r(l:2:zero_fill(l)*2,l); 
tphasel(:,l)  =  phase_r(l:2:zero_fill(l)*2,l); 
for  col  =  2:size(X_r,2) 

tmagcol(:,col-l)  =  mag_r(l:2:zero_fill(col)*2,col); 
tphasecol(:, col-1)  =  phase_r(l:2:zero_fill(col)*2,col); 
end 

mag_r_vecl(:,l)  =  [tmagl(:,l);  mag_r(zero_fill(l)*2+l:length(mag_r),l)]; 
phase_r_vecl(:,l)  =  [tphasel(:,l);  phase_r(zero_fill(l)*2+l:length(phase_r),l)]; 
for  col  =  2:size(X_r,2) 
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[tmagcol(:,col-l); 


mag_matcol(:, col-1)  = 

mag_r(zero_fill(col)*2+l:length(mag_r),col)]; 

phase_matcol(:, col-1)  =  [tphasecol(:,col-l); 

phase_r(zero_fill(col)*2+ 1  :length(phase_r),col)]; 

end 

magrvcccol  =  reshape(mag_matcol,  size(mag_matcol,l)*size(mag_matcol,2), 

D; 

phaserveccol  =  reshape(phase_matcol, 

size(phase_matcol,  l)*size(phase_matcol,2),  1); 

magrvccshort  =  [mag_r_vecl;  mag_r_veccol]; 

phaservecshort  =  [phaservec 1 ;  phase  r  veccol]; 


[inphase  r,  quadr]  =  pol2cart(phase_r_vec_short,  mag  r  vec  short); 
%converts  the  mag  and  phase  values  to  in  phase  and  quadriture  components 
%for  QAM  demodulation 

minmax  =  [min(inphaser),  max(inphaser);  min(quad  r),  max(quadr)]; 
code  =  qaskdecomod(inphase_r,  quad  r,  Q,  minmax); 

%decodes  the  QAM  symbols  into  integers  0:Q-1 


a.  qaskdecomod.m 

function  [msg]=qaskdecomod(x,y,m, minmax) 

%QASKDECO  Decodes  QASK  mapped  signal  back  to  integer  symbols. 

%  MSG  =  QASKDECO(X,  Y,  M)  decodes  the  message  signal  MSG 
from  the 

%  in-phase  component  X  and  quadrature  component  Y  with  the  M-ary 

%  number  M.  M  must  equal  2AK  with  K  being  an  positive  integer.  The 

%  minimum  and  the  maximum  numbers  are  assumed  to  be  the  same  as 
in  the 

%  output  from  QASKENCO. 

% 

%  MSG  =  QASKDECO(X,  Y,  M,  MINMAX)  decodes  the  information 
where 

%  the  maximum  and  minimum  values  of  X  and  Y  are  given  in 

MINMAX 
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% 

by  the  form: 

% 

X  min 

X  max 

% 

MINMAX  =  | 

% 

Y  min 

Y  max 

% 

% 

See  also  QASKENCO. 

% 

Last  modified  5/11/01 

scale  =  0;  %logic  statement  to  scale  received  symbols,  1  =  yes,  o/w  = 
no 

%if  scale  ==  0 

%  disp(’NO  SCALING  OF  RECEIVED  QAM  CONSTELLATION 
PERFORMED  IN  QAM  DECODER’) 

%end 

error(nargchk(  1 ,4,nargin)); 
if  (nargin  <=  2)  |  (nargin  >  4) 

error('Wrong  number  of  input  arguments.’) 
end; 

M  =  m; 
if  isstr(m) 

error) [’QASK  constellation  has  changed  to  use  Gray  code. ',  ... 
sprintf('\n’),... 

'Use  COMMUPDT ("FILENAME")  to  update  your  SIMULINK 
model’... 

sprintf('\n’),... 

’before  the  simulation.’]); 

end; 

K  =  log(m)  ./  log(2); 
if  m  ~=  2AK 

error(’M  must  equal  to  2AK  where  K  is  an  integer.') 
end; 

xx  =  constlay(K,  1); 

[leny,  lenx]  =  size(xx); 
if  (nargin  ==  4) 
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[xmat,ymat]  =  qaskenco(M); 

%xmat  and  ymat  are  the  inphase  and  quad  components,  respectively,  of 
the  Mary  QAM 

%constellation 

MAT  =  [xmat.';ymat.']; 

%2  by  M  matrix  of  vectors  in  2D  representing  the  M  Qam  symbols  in 
the  constelation 

for  k=l:length(x) 

E  =  [x(k);  y(k)]*ones(l,length(xmat))  -  MAT; 

[v,  index]  =  min([l,l]*abs(E)); 
msg(k)  =  index- 1 ; 
end 
end; 

[x_m,  x_n]  =  size(x); 

[yjn,  y_n]  =  size(y); 
if  min(x_m,  x_n)  *  min(y_m,  y_n)  ~=  1 
error(’Input  X  and  Y  must  be  vectors.'); 
end; 

if  x_m  >  x_n 
msg  =  msg'; 
end; 

7.  deinterleaver.m 

function  deint_mat  =  deinterleaver(code,  N) 

%Last  modified  3/19/01 
rows  =  length(code)/N; 

%number  of  rows  in  matrix  for  RS  decoding 
deint_mat  =  reshape(code,  rows,  N); 

%matrix  of  symbols  deinterleaved,  ready  for  RS  decoder 
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8.  rsdecoder.m 


function  decode  =  rsdecoder(code_deint,  N,  k,  zero_symbols) 

%Last  modified  3/19/01 

deco  =  rsdeco(code_deint,  N,  k,  'rs/decimal'); 

%message  length  k  columns  matrix 

deco2  =  reshape(deco,  1,  size(deco,l)*size(deco,2)); 

%converts  to  row  vector 

decode  =  deco2(zero_symbols+l:length(deco2)); 

%removes  zeroes  added  in  transmitter 

%zero  symbols  passed  from  xmitter  function  reedsolomon.m 


9.  error  check.m 

function  ber  =  error_check(bs,  bsr) 

%This  function  determines  the  Bit  Error  Rate  (BER)  for  the  system. 
%last  modified  7/2/01 
be  =  0; 

for  1  =  1  dength(bs) 
be  =  be  +  (bs(l)~=bsr(l)); 
end 

ber  =  be/1; 


C.  PERFORMANCE  ANALYSIS  CODE 

1.  ofdmsimcontrolfunc.m 

function  [BER,  SNR  est]  =  ofdm_sim_control_func(SNR,  fc,  R,  rough,  range, 
depth) 

%Control  m- file  for  ofdm  simulation 
fn  =  0; 

modulate  =  2;  %DSBSC  modulation  logic  control  value 

noise  =  1 ; 
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for  avn  =1:10 
count  =  0; 

for  noisefactor  =  .15:. 01:. 65 
count  =  count  +  1 ; 

synchdec  =  0;  %synchronization  logic  control  value,  0  ==  NO,  1 

==  YES 

[zero_fill,  Q,  N,  data_carriers,  N_bar,  N_subc_act,  zero_symbols,  k,  bits, 
info_sig,  code_stream,  bit_stream,  X,  mod_xmit_sig,  ... 

Fs_prime,  Fs  xmit,  fn,  blocklength,  zerofraction,  SNR  awgn,  chanapp, 
W,  n_lpf2,  lenxint,  F,  q,  codebook]  =  ... 

ofdm_sim_xmitter_fimc(synchdec,  noise,  noise_factor,  fn,  R,  fc,  SNR); 

[SNR_est(count,avn),  BER(count,avn)]  =  ofdm_sim_receiver_func(W, 
zero  fill,  Q,  N,  data_carriers,  N_bar,  N_subc_act,  zero_symbols,  k,  bits,  info_sig, 

code_stream,  bit_stream,  X,  mod_xmit_sig,  Fs_prime,  Fs_xmit,  fn, 
block  length,  zero  fraction,  SNR  awgn,  ... 

depth,  range,  rough,  chan  app,  modulate,  noise,  fc,  synchdec,  n_lpf2, 
len  xint,  F,  q,  codebook); 

end 

end 


2.  ofdmperformance.m 

%OFDM  performance  analysis  m-file 

%runs  ofdm  sim  control  func  m-file  over  all  mmpe  channels 

%last  modified  6/28/01 


warning  off 
clear  all 


%%% 


%8-ary  QAM  PERFORMANCE 

SNR  =  17 

fc  =  8e3 

R=5e3 

count  =  0; 
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for  range  =  0:1  %0  ==  2  km,  1  ==  4  km 

count  =  count  +  1 ; 
rough  =  0; 
depth  =  0; 

[BERw,  SNR_estw]  =  ofdm_sim_control_func(SNR,  fc,  R,  rough,  range, 
depth); 

BER(:,:, count)  =  BERw; 

SNR_est(:,:, count)  =  SNR_estw; 
end 

clear  BERw  SNR_estw 
fc  =  10e3 
R  =  6e3 

for  range  =  0:1  %0  ==  2  km,  1  ==  4  km 

count  =  count  +  1 ; 
rough  =  0; 
depth  =  0; 

[BERw,  SNR_estw]  =  ofdm_sim_control_func(SNR,  fc,  R,  rough,  range, 
depth); 

BER(:,:, count)  =  BERw; 

SNR_est(:,:, count)  =  SNR_estw; 
end 

clear  BERw  SNR_estw 
fc  =  12e3 
R=  7.5e3 

for  range  =  0:1  %0  ==  2  km,  1  ==  4  km 

count  =  count  +  1 ; 
rough  =  0; 
depth  =  0; 

[BERw,  SNR_estw]  =  ofdm_sim_control_func(SNR,  fc,  R,  rough,  range, 
depth); 

BER(:,:, count)  =  BERw; 

SNR_est(:,:, count)  =  SNR_estw; 
end 
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clear  BERw  SNR  estw 


%%% 

%  16-ary  QAM  PERFORMANCE 
SNR  =  21 
fc  =  6e3 
R=  5e3 
depth  =  0; 

for  range  =  0:1  %0  ==  2  km,  1  ==  4  km 

for  rough  =  0:2:4 
count  =  count  +  1 ; 

[BERw,  SNR_estw]  =  ofdm_sim_control_func(SNR,  fc,  R,  rough,  range, 
depth); 

BER(:,:, count)  =  BERw; 

SNR_est(:,:, count)  =  SNR_estw; 
end  %rough 
end  %range 
clear  BERw  SNR_estw 
depth  =  1 ; 

for  range  =  0: 1  %0  ==  2  km,  1  ==  4  km 

for  rough  =  0:2:2 
count  =  count  +  1 ; 

[BERw,  SNR_estw]  =  ofdm_sim_control_func(SNR,  fc,  R,  rough,  range, 
depth); 

BER(:,:, count)  =  BERw; 

SNR_est(:,:, count)  =  SNR_estw; 
end  %rough 
end  %range 
clear  BERw  SNR_estw 
fc  =  12e3 
R  =  10e3 

for  range  =  0:1  %0  ==  2  km,  1  ==  4  km 

count  =  count  +  1 ; 
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rough  =  0; 
depth  =  0; 

[BERw,  SNR_estw]  =  ofdm_sim_control_func(SNR,  fc,  R,  rough,  range, 
depth); 

BER(:,:, count)  =  BERw; 

SNR_est(:,:, count)  =  SNR_estw; 
end 

clear  BERw  SNR_estw 
fc  =  8e3 

R  =  (6+2/3)*  le3 

for  range  =  0:1  %0  ==  2  km,  1  ==  4  km 

count  =  count  +  1 ; 
rough  =  0; 
depth  =  0; 

[BERw,  SNR_estw]  =  ofdm_sim_control_func(SNR,  fc,  R,  rough,  range, 
depth); 

BER(:,:, count)  =  BERw; 

SNR_est(:,:, count)  =  SNR_estw; 
end 

clear  BERw  SNR_estw 
fc  =  10e3 
R=  (8+0/3)*  le3 

for  range  =  0:1  %0  ==  2  km,  1  ==  4  km 

count  =  count  +  1 ; 
rough  =  0; 
depth  =  0; 

[BERw,  SNR_estw]  =  ofdm_sim_control_func(SNR,  fc,  R,  rough,  range, 
depth); 

BER(:,:, count)  =  BERw; 

SNR_est(:,:, count)  =  SNR_estw; 
end 

clear  BERw  SNR  estw 
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%%% 

%32-ary  QAM  PERFORMANCE 
SNR  =  25 
fc  =  10e3 
R  =  (10+0/3)*  le3 

for  range  =  0:1  %0  ==  2  km,  1  ==  4  km 

count  =  count  +  1 ; 
rough  =  0; 
depth  =  0; 

[BERw,  SNR_estw]  =  ofdm_sim_control_func(SNR,  fc,  R,  rough,  range, 
depth); 

BER(:,:, count)  =  BERw; 

SNR_est(:,:, count)  =  SNR_estw; 
end 

clear  BERw  SNR_estw 
fc  =  8e3 

R  =  (8+l/3)*le3 

for  range  =  0:1  %0  ==  2  km,  1  ==  4  km 

count  =  count  +  1 ; 
rough  =  0; 
depth  =  0; 

[BERw,  SNR_estw]  =  ofdm_sim_control_func(SNR,  fc,  R,  rough,  range, 
depth); 

BER(:,:, count)  =  BERw; 

SNR_est(:,:, count)  =  SNR_estw; 
end 

clear  BERw  SNR_estw 
fc  =  6e3 

R=  (5+0/3)*  le3 
depth  =  0; 

for  range  =  0:1  %0  ==  2  km,  1  ==  4  km 

for  rough  =  0:2:4 
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count  =  count  +  1 ; 

[BERw,  SNR_estw]  =  ofdm_sim_control_func(SNR,  fc,  R,  rough,  range, 
depth); 

BER(:,:, count)  =  BERw; 

SNR_est(:,:, count)  =  SNR_estw; 
end  %rough 
end  %range 
clear  BERw  SNR_estw 
fc  =  6e3 

R  =  (5+0/3)*  le3 
depth  =  1 ; 

for  range  =  0:1  %0  ==  2  km,  1  ==  4  km 

for  rough  =  0:2:2 
count  =  count  +  1 ; 

[BERw,  SNR_estw]  =  ofdm_sim_control_func(SNR,  fc,  R,  rough,  range, 
depth); 

BER(:,:, count)  =  BERw; 

SNR_est(:,:, count)  =  SNR_estw; 
end  %rough 
end  %range 
clear  BERw  SNR_estw 

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 


%64ary  QAM  PERFORMANCE 
SNR  =  27 
fc  =  8e3 


R=  (10+0/3)*  le3 

for  range  =  0:1  %0  ==  2  km,  1  ==  4  km 

count  =  count  +  1 ; 
rough  =  0; 
depth  =  0; 

[BERw,  SNR_estw]  =  ofdm_sim_control_func(SNR,  fc,  R,  rough,  range, 
depth); 
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BER(:,:, count)  =  BERw; 

SNR_est(:,:, count)  =  SNR_estw; 
end 

clear  BERw  SNR_estw 
fc  =  6e3 

R  =  (5+0/3)*  le3 
depth  =  0; 

for  range  =  0:1  %0  ==  2  km,  1  ==  4  km 

for  rough  =  0:2:4 
count  =  count  +  1 ; 

[BERw,  SNR_estw]  =  ofdm_sim_control_func(SNR,  fc,  R,  rough,  range, 
depth); 

BER(:,:, count)  =  BERw; 

SNR_est(:,:, count)  =  SNR_estw; 
end  %rough 
end  %range 
clear  BERw  SNR_estw 
fc  =  6e3 

R  =  (5+0/3)*  1  e3 
depth  =  1 ; 

for  range  =  0:1  %0  ==  2  km,  1  ==  4  km 

for  rough  =  0:2:2 
count  =  count  +  1 ; 

[BERw,  SNR_estw]  =  ofdm_sim_control_func(SNR,  fc,  R,  rough,  range, 
depth); 

BER(:,:, count)  =  BERw; 

SNR_est(:,:, count)  =  SNR_estw; 
end  %rough 
end  %range 
clear  BERw  SNR_estw 
warning  on 

save  perform  BER  SNR_est 
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3.  performplotter.m 


%m-file  to  plot  BER  vs  SNR_est  for  performance  results 
clear  all 
load  perform 
fn  =  0; 


close  all 


%%%%%%%%%%%%%%%%%%%%%%%%%%% 


%8  QAM 
for  1  =  1 :6 

mn_ber8(:,l)  =  mean(BER(:,:,l),2); 
mn_snr8(:,l)  =  mean(SNR_est(:,:,l),2); 
end 
fnc  =  0 
fnc  =  fnc+ 1 
figure(fnc) 

plot(mn_snr8,mn_ber8) 
title(’8  QAM  Composite  Results') 
xlabel('SNR  (dB)’),  ylabel(’BER’) 
legend(T,  '2',  '3',  '4',  '5',  '6') 

%%%%%%%%%%%%%%%%%%%%%%%%% 


%16  QAM 


for  1  =  7:22 

mn_berl6(:,l)  =  mean(BER(:,:,l),2); 
mn_snrl6(:,l)  =  mean(SNR_est(:,:,l),2); 
end 

fnc  =  fnc+ 1 
figure(fnc) 

plot(mn_snrl6,mn_berl6) 
title(’16  QAM  Composite  Results') 
xlabel('SNR  (dB)’),  ylabel(’BER’) 


legend(T,  '2',  '3',  '4',  ’5’,  '6',  7’,  '8',  '9',  '10',  '1 1’,  '12',  '13',  '14',  '15',  '16') 


%%%%%%%%%%%%%%%%%%%%%%%%%%% 
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%32  QAM 


for  1  =  23:36 

mn_ber32(:,l)  =  mean(BER(:,:,l),2); 
mn_snr32(:,l)  =  mean(SNR_est(:,:,l),2); 
end 

fnc  =  fnc+ 1 
figure(fnc) 

plot(mn_snr32,mn_ber32) 
title(’32  QAM  Composite  Results') 
xlabel('SNR  (dB)’),  ylabel(’BER’) 


legend(T,  '2',  '3',  '4',  '5',  '6',  7’,  '8',  '9',  '10',  '1 1’,  '12',  '13',  '14') 


%%%%%%%%%%%%%%%%%%%%%%%%% 


%64  QAM 


for  1  =  37:48 

mn_ber64(:,l)  =  mean(BER(:,:,l),2); 
mn_snr64(:,l)  =  mean(SNR_est(:,:,l),2); 
end 

fnc  =  fnc+ 1 
figure(fnc) 

plot(mn_snr64,mn_ber64) 
title('64  QAM  Composite  Results') 
xlabel('SNR  (dB)’),  ylabel(’BER’) 


legend(T,  '2',  '3',  '4',  '5',  '6',  7’,  '8',  '9',  '10',  '1 1’ 


'12') 


%break 


%%% 


%8-ary  QAM  PERFORMANCE 
SNR  =  17 
fc  =  8e3 
R=5e3 


%%% 

fn  =  fin  +  1 ;  figure(fn),  orient  tall 
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plot(SNR_est(:,:,fn),  BER(:,:,fn),  'b.’) 
hold  on 

plot(mean(SNR_est(:,:,fn),2),  mean(BER(:,:,fn),2),  'r-') 
hold  off 

title(['SNR  =  ’,num2str(SNR),  'dB,  f_c  =  ',num2str(fc/1000),'kHz,  R 

',num2str(R/l 000), ’kbps,  Range  =  2  km,  Rough  =  0,  Water  Depth  =  100  m']) 

xlabel(’SNR  (dB)'),  ylabel(’BER’) 

%%% 

fn  =  fin  +  1 ;  figure(fn),  orient  tall 
plot(SNR_est(:,:,fn),  BER(:,:,fn),  'b.') 
hold  on 

plot(mean(SNR_est(:,:,fn),2),  mean(BER(:,:,fn),2),  ’r-’) 
hold  off 

title(['SNR  =  ’,num2str(SNR),  ’dB,  f_c  =  ',num2str(fc/1000),'kHz,  R 

’,num2str(R/l  000), 'kbps,  Range  =  4  km,  Rough  =  0,  Water  Depth  =  100  m']) 

xlabel('SNR  (dB)’),  ylabel(’BER’) 

%%% 

fc  =  10e3 

R  =  6e3 

%%% 

fn  =  fn  +  1 ;  figure(fn),  orient  tall 
plot(SNR_est(:,:,fn),  BER(:,:,fn),  ’b.’) 
hold  on 

plot(mean(SNR_est(:,:,fn),2),  mean(BER(:,:,fn),2),  ’r-’) 
hold  off 

title(['SNR  =  ’,num2str(SNR),  'dB,  f_c  =  ',num2str(fc/1000),'kHz,  R 
’,num2str(R/l  000), 'kbps,  Range  =  2  km,  Rough  =  0,  Water  Depth  =  100  m']) 

xlabel('SNR  (dB)’),  ylabel(’BER’) 

%%% 

fn  =  fn  +  1 ;  figure(fn),  orient  tall 
plot(SNR_est(:,:,fn),  BER(:,:,fn),  'b.') 
hold  on 

plot(mean(SNR_est(:,:,fn),2),  mean(BER(:,:,fn),2),  'r-') 
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hold  off 

title(['SNR  =  ’,num2str(SNR),  ’dB,  f_c  =  ’,num2str(fc/1000),'kHz,  R  = 
',num2str(R/l  000), 'kbps,  Range  =  4  km,  Rough  =  0,  Water  Depth  =  100  m']) 

xlabel('SNR  (dB)’),  ylabel(’BER’) 

%%% 

fc  =  12e3 

R  =  7.5e3 

%%% 

fn  =  fn  +  1 ;  figure(fn),  orient  tall 
plot(SNR_est(:,:,fn),  BER(:,:,fn),  'b.') 
hold  on 

plot(mean(SNR_est(:,:,fn),2),  mean(BER(:,:,fn),2),  ’r-’) 
hold  off 

title(['SNR  =  ’,num2str(SNR),  'dB,  f_c  =  ',num2str(fc/1000),'kHz,  R  = 

’,num2str(R/l  000), 'kbps,  Range  =  2  km,  Rough  =  0,  Water  Depth  =  100  m']) 

xlabel('SNR  (dB)’),  ylabel(’BER’) 

%%% 

fn  =  fn  +  1 ;  figure(fn),  orient  tall 
plot(SNR_est(:,:,fn),  BER(:,:,fn),  'b.') 
hold  on 

plot(mean(SNR_est(:,:,fn),2),  mean(BER(:,:,fn),2),  'r-') 
hold  off 

title(['SNR  =  ’,num2str(SNR),  'dB,  f_c  =  ',num2str(fc/1000),'kHz,  R  = 

’,num2str(R/l  000), 'kbps,  Range  =  4  km,  Rough  =  0,  Water  Depth  =  100  m']) 

xlabel('SNR  (dB)’),  ylabel(’BER’) 

%%% 


%%% 


%  16-ary  QAM  PERFORMANCE 
SNR  =  21 
fc  =  6e3 
R=  5e3 
%%% 
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fn  =  fn  +  1 ;  figure(fn),  orient  tall 
plot(SNR_est(:,:,fn),  BER(:,:,fn),  ’b.’) 
hold  on 

plot(mean(SNR_est(:,:,fn),2),  mean(BER(:,:,fn),2),  'r-') 
hold  off 

title(['SNR  =  ’,num2str(SNR),  'dB,  f_c  =  ',num2str(fc/1000),'kHz,  R 
',num2str(R/l  000), 'kbps,  Range  =  2  km,  Rough  =  0,  Water  Depth  =  100  m']) 

xlabel('SNR  (dB)’),  ylabel(’BER’) 

%%% 

fn  =  fn  +  1 ;  figure(fn),  orient  tall 
plot(SNR_est(:,:,fn),  BER(:,:,fn),  'b.') 
hold  on 

plot(mean(SNR_est(:,:,fn),2),  mean(BER(:,:,fn),2),  'r-') 
hold  off 

title(['SNR  =  ’,num2str(SNR),  'dB,  f_c  =  ',num2str(fc/1000),'kHz,  R 

',num2str(R/l 000), ’kbps,  Range  =  2  km,  Rough  =  2,  Water  Depth  =  100  m']) 

xlabel(’SNR  (dB)'),  ylabel(’BER’) 

%%% 

fn  =  fn  +  1 ;  figure(fn),  orient  tall 
plot(SNR_est(:,:,fn),  BER(:,:,fn),  'b.') 
hold  on 

plot(mean(SNR_est(:,:,fn),2),  mean(BER(:,:,fn),2),  ’r-’) 
hold  off 

title(['SNR  =  ’,num2str(SNR),  ’dB,  f_c  =  ',num2str(fc/1000),'kHz,  R 

’,num2str(R/l  000), 'kbps,  Range  =  2  km,  Rough  =  4,  Water  Depth  =  100  m']) 

xlabel('SNR  (dB)’),  ylabel(’BER’) 

%%% 

fn  =  fn  +  1 ;  figure(fn),  orient  tall 
plot(SNR_est(:,:,fn),  BER(:,:,fn),  'b.') 
hold  on 

plot(mean(SNR_est(:,:,fn),2),  mean(BER(:,:,fn),2),  ’r-’) 
hold  off 

title(['SNR  =  ’,num2str(SNR),  'dB,  f_c  =  ',num2str(fc/1000),'kHz,  R 

’,num2str(R/l  000), 'kbps,  Range  =  4  km,  Rough  =  0,  Water  Depth  =  100  m']) 
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xlabel('SNR  (dB)’),  ylabel(’BER’) 

%%% 

fn  =  fn  +  1 ;  figure(fn),  orient  tall 
plot(SNR_est(:,:,fn),  BER(:,:,fn),  'b.’) 
hold  on 

plot(mean(SNR_est(:,:,fn),2),  mean(BER(:,:,fn),2),  'r-') 
hold  off 

title(['SNR  =  ’,num2str(SNR),  ’dB,  f_c  =  ’,num2str(fc/1000),’kHz,  R  = 

’,num2str(R/l  000), 'kbps,  Range  =  4  km,  Rough  =  2,  Water  Depth  =  100  m']) 

xlabel('SNR  (dB)’),  ylabel(’BER’) 

%%% 

fn  =  fn  +  1 ;  figure(fn),  orient  tall 
plot(SNR_est(:,:,fn),  BER(:,:,fn),  ’b.’) 
hold  on 

plot(mean(SNR_est(:,:,fn),2),  mean(BER(:,:,fn),2),  'r-') 
hold  off 

title(['SNR  =  ’,num2str(SNR),  'dB,  f_c  =  ’,num2str(fc/1000),'kHz,  R  = 

’,num2str(R/l  000), 'kbps,  Range  =  4  km,  Rough  =  4,  Water  Depth  =  100  m']) 

xlabel('SNR  (dB)’),  ylabel(’BER’) 

%%% 

%depth  =  0; 

%for  range  =  0:1  %0  ==  2  km,  1  ==  4  km 

%  for  rough  =  0:2:4 
%  count  =  count  +  1 ; 

%  [BERw,  SNR_estw]  =  ofdm_sim_control_func(SNR,  fc,  R,  rough,  range, 
depth); 

%  BER(:,:, count)  =  BERw; 

%  SNR_est(:,:, count)  =  SNR_estw; 

%end  %rough 
%end  %range 
%clear  BERw  SNR_estw 
%%% 

fn  =  fn  +  1 ;  figure(fn),  orient  tall 
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plot(SNR_est(:,:,fn),  BER(:,:,fn),  'b.’) 
hold  on 

plot(mean(SNR_est(:,:,fn),2),  mean(BER(:,:,fn),2),  'r-') 
hold  off 

title(['SNR  =  ’,num2str(SNR),  'dB,  f_c  =  ',num2str(fc/1000),’kHz,  R 
',num2str(R/l 000), ’kbps,  Range  =  2  km,  Rough  =  0,  Water  Depth  =  340  m']) 

xlabel(’SNR  (dB)'),  ylabel(’BER’) 

%%% 

fn  =  fin  +  1 ;  figure(fn),  orient  tall 
plot(SNR_est(:,:,fn),  BER(:,:,fn),  'b.') 
hold  on 

plot(mean(SNR_est(:,:,fn),2),  mean(BER(:,:,fn),2),  ’r-’) 
hold  off 

title(['SNR  =  ’,num2str(SNR),  ’dB,  f_c  =  ',num2str(fc/1000),'kHz,  R 

’,num2str(R/l  000), 'kbps,  Range  =  2  km,  Rough  =  2,  Water  Depth  =  340  m']) 

xlabel('SNR  (dB)’),  ylabel(’BER’) 

%%% 

fn  =  fn  +  1 ;  figure(fn),  orient  tall 
plot(SNR_est(:,:,fn),  BER(:,:,fn),  ’b.’) 
hold  on 

plot(mean(SNR_est(:,:,fn),2),  mean(BER(:,:,fn),2),  ’r-’) 
hold  off 

title(['SNR  =  ’,num2str(SNR),  ’dB,  f_c  =  ',num2str(fc/1000),'kHz,  R 

’,num2str(R/l  000), 'kbps,  Range  =  4  km,  Rough  =  0,  Water  Depth  =  340  m']) 

xlabel('SNR  (dB)’),  ylabel(’BER’) 

%%% 

fn  =  fn  +  1 ;  figure(fn),  orient  tall 
plot(SNR_est(:,:,fn),  BER(:,:,fn),  'b.') 
hold  on 

plot(mean(SNR_est(:,:,fn),2),  mean(BER(:,:,fn),2),  'r-') 
hold  off 

title(['SNR  =  ’,num2str(SNR),  ’dB,  f_c  =  ’,num2str(fc/1000),'kHz,  R 

',num2str(R/l 000), 'kbps,  Range  =  4  km,  Rough  =  2,  Water  Depth  =  340  m’]) 

xlabel('SNR  (dB)’),  ylabel(’BER’) 
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%%% 
fc  =  12e3 
R  =  10e3 

%%% 

fn  =  fn  +  1 ;  figure(fn),  orient  tall 
plot(SNR_est(:,:,fn),  BER(:,:,fn),  ’b.') 
hold  on 

plot(mean(SNR_est(:,:,fn),2),  mean(BER(:,:,fn),2),  ’r-’) 
hold  off 

title(['SNR  =  ’,num2str(SNR),  ’dB,  f_c  =  ',num2str(fc/1000),'kHz,  R 

',num2str(R/l  000), 'kbps,  Range  =  2  km,  Rough  =  0,  Water  Depth  =  100  m']) 

xlabel('SNR  (dB)’),  ylabel(’BER’) 

%%% 

fn  =  fn  +  1 ;  figure(fn),  orient  tall 
plot(SNR_est(:,:,fn),  BER(:,:,fn),  'b.') 
hold  on 

plot(mean(SNR_est(:,:,fn),2),  mean(BER(:,:,fn),2),  ’r-’) 
hold  off 

title(['SNR  =  ’,num2str(SNR),  'dB,  f_c  =  ',num2str(fc/1000),'kHz,  R 

',num2str(R/l  000), 'kbps,  Range  =  4  km,  Rough  =  0,  Water  Depth  =  100  m']) 

xlabel('SNR  (dB)’),  ylabel(’BER’) 

%%% 

fc  =  8e3 

R  =  (6+2/3)*  le3 

%%% 

fn  =  fn  +  1 ;  figure(fn),  orient  tall 
plot(SNR_est(:,:,fn),  BER(:,:,fn),  'b.') 
hold  on 

plot(mean(SNR_est(:,:,fn),2),  mean(BER(:,:,fn),2),  'r-') 
hold  off 

title(['SNR  =  ’,num2str(SNR),  'dB,  f_c  =  ',num2str(fc/1000),'kHz,  R 
’,num2str(R/l  000), 'kbps,  Range  =  2  km,  Rough  =  0,  Water  Depth  =  100  m']) 

xlabel('SNR  (dB)’),  ylabel(’BER’) 
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%%% 

fn  =  fn  +  1 ;  figure(fn),  orient  tall 
plot(SNR_est(:,:,fn),  BER(:,:,fn),  ’b.’) 
hold  on 

plot(mean(SNR_est(:,:,fn),2),  mean(BER(:,:,fn),2),  'r-') 
hold  off 

title(['SNR  =  ’,num2str(SNR),  ’dB,  f_c  =  ’,num2str(fc/1000),'kHz,  R  = 
’,num2str(R/l  000), 'kbps,  Range  =  4  km,  Rough  =  0,  Water  Depth  =  100  m']) 

xlabel('SNR  (dB)’),  ylabel(’BER’) 

%%% 

fc  =  10e3 

R  =  (8+0/3)*  1  e3 

%%% 

fn  =  fn  +  1 ;  figure(fn),  orient  tall 
plot(SNR_est(:,:,fn),  BER(:,:,fn),  'b.') 
hold  on 

plot(mean(SNR_est(:,:,fn),2),  mean(BER(:,:,fn),2),  ’r-’) 
hold  off 

title(['SNR  =  ’,num2str(SNR),  'dB,  f_c  =  ',num2str(fc/1000),'kHz,  R  = 

',num2str(R/l  000), 'kbps,  Range  =  2  km,  Rough  =  0,  Water  Depth  =  100  m']) 

xlabel('SNR  (dB)’),  ylabel(’BER’) 

%%% 

fn  =  fn  +  1 ;  figure(fn),  orient  tall 
plot(SNR_est(:,:,fn),  BER(:,:,fn),  'b.') 
hold  on 

plot(mean(SNR_est(:,:,fn),2),  mean(BER(:,:,fn),2),  'r-') 
hold  off 

title(['SNR  =  ’,num2str(SNR),  'dB,  f_c  =  ’,num2str(fc/1000),'kHz,  R  = 

’,num2str(R/l  000), 'kbps,  Range  =  4  km,  Rough  =  0,  Water  Depth  =  100  m']) 

xlabel('SNR  (dB)’),  ylabel(’BER’) 

%%% 

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 
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%32-ary  QAM  PERFORMANCE 
SNR  =  25 
fc  =  10e3 
R=  (10+0/3)*  le3 

%%% 

fn  =  fn  +  1 ;  figure(fn),  orient  tall 
plot(SNR_est(:,:,fn),  BER(:,:,fn),  ’b.’) 
hold  on 

plot(mean(SNR_est(:,:,fn),2),  mean(BER(:,:,fn),2),  'r-') 
hold  off 

title(['SNR  =  ’,num2str(SNR),  ’dB,  f_c  =  ',num2str(fc/1000),'kHz,  R 

’,num2str(R/l  000), 'kbps,  Range  =  2  km,  Rough  =  0,  Water  Depth  =  100  m']) 

xlabel('SNR  (dB)’),  ylabel(’BER’) 

%%% 

fn  =  fn  +  1 ;  figure(fn),  orient  tall 
plot(SNR_est(:,:,fn),  BER(:,:,fn),  'b.') 
hold  on 

plot(mean(SNR_est(:,:,fn),2),  mean(BER(:,:,fn),2),  'r-') 
hold  off 

title(['SNR  =  ’,num2str(SNR),  'dB,  f_c  =  ',num2str(fc/1000),'kHz,  R 

',num2str(R/l 000), 'kbps,  Range  =  4  km,  Rough  =  0,  Water  Depth  =  100  m']) 

xlabel('SNR  (dB)’),  ylabel(’BER’) 

%%% 

fc  =  8e3 

R  =  (8+l/3)*le3 
%%% 

fn  =  fn  +  1 ;  figure(fn),  orient  tall 
plot(SNR_est(:,:,fn),  BER(:,:,fn),  'b.') 
hold  on 

plot(mean(SNR_est(:,:,fn),2),  mean(BER(:,:,fn),2),  'r-') 
hold  off 

title(['SNR  =  ’,num2str(SNR),  'dB,  f_c  =  ',num2str(fc/1000),'kHz,  R 
’,num2str(R/l  000), 'kbps,  Range  =  2  km,  Rough  =  0,  Water  Depth  =  100  m']) 
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xlabel('SNR  (dB)’),  ylabel(’BER’) 

%%% 

fn  =  fn  +  1 ;  figure(fn),  orient  tall 
plot(SNR_est(:,:,fn),  BER(:,:,fn),  'b.’) 
hold  on 

plot(mean(SNR_est(:,:,fn),2),  mean(BER(:,:,fn),2),  'r-') 
hold  off 

title(['SNR  =  ’,num2str(SNR),  ’dB,  f_c  =  ’,num2str(fc/1000),’kHz,  R 
’,num2str(R/l  000), 'kbps,  Range  =  4  km,  Rough  =  0,  Water  Depth  =  100  m']) 

xlabel('SNR  (dB)'),  ylabel(’BER’) 

%%% 

fc  =  6e3 

R  =  (5+0/3)*  le3 

%%% 

fn  =  fn  +  1 ;  figure(fn),  orient  tall 
plot(SNR_est(:,:,fn),  BER(:,:,fn),  'b.') 
hold  on 

plot(mean(SNR_est(:,:,fn),2),  mean(BER(:,:,fn),2),  'r-') 
hold  off 

title(['SNR  =  ’,num2str(SNR),  'dB,  f_c  =  ',num2str(fc/1000),'kHz,  R 

',num2str(R/l 000), 'kbps,  Range  =  2  km,  Rough  =  0,  Water  Depth  =  100  m']) 

xlabel('SNR  (dB)’),  ylabel(’BER’) 

%%% 

fn  =  fn  +  1 ;  figure(fn),  orient  tall 
plot(SNR_est(:,:,fn),  BER(:,:,fn),  'b.') 
hold  on 

plot(mean(SNR_est(:,:,fn),2),  mean(BER(:,:,fn),2),  'r-') 
hold  off 

title(['SNR  =  ’,num2str(SNR),  'dB,  f_c  =  ',num2str(fc/1000),'kHz,  R 

',num2str(R/l 000), 'kbps,  Range  =  2  km,  Rough  =  2,  Water  Depth  =  100  m']) 

xlabel('SNR  (dB)’),  ylabel(’BER’) 

%%% 

fn  =  fn  +  1 ;  figure(fn),  orient  tall 
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plot(SNR_est(:,:,fn),  BER(:,:,fn),  'b.’) 
hold  on 

plot(mean(SNR_est(:,:,fn),2),  mean(BER(:,:,fn),2),  'r-') 
hold  off 

title(['SNR  =  ’,num2str(SNR),  'dB,  f_c  =  ’,num2str(fc/1000),'kHz,  R 
',num2str(R/l 000), ’kbps,  Range  =  2  km,  Rough  =  4,  Water  Depth  =  100  m']) 

xlabel(’SNR  (dB)'),  ylabel(’BER’) 

%%% 

fn  =  fin  +  1 ;  figure(fn),  orient  tall 
plot(SNR_est(:,:,fn),  BER(:,:,fn),  'b.') 
hold  on 

plot(mean(SNR_est(:,:,fn),2),  mean(BER(:,:,fn),2),  ’r-’) 
hold  off 

title(['SNR  =  ’,num2str(SNR),  ’dB,  f_c  =  ',num2str(fc/1000),'kHz,  R 

’,num2str(R/l  000), 'kbps,  Range  =  4  km,  Rough  =  0,  Water  Depth  =  100  m']) 

xlabel('SNR  (dB)’),  ylabel(’BER’) 

%%% 

fn  =  fn  +  1 ;  figure(fn),  orient  tall 
plot(SNR_est(:,:,fn),  BER(:,:,fn),  ’b.’) 
hold  on 

plot(mean(SNR_est(:,:,fn),2),  mean(BER(:,:,fn),2),  ’r-’) 
hold  off 

title(['SNR  =  ’,num2str(SNR),  ’dB,  f_c  =  ',num2str(fc/1000),'kHz,  R 

’,num2str(R/l  000), 'kbps,  Range  =  4  km,  Rough  =  2,  Water  Depth  =  100  m']) 

xlabel('SNR  (dB)’),  ylabel(’BER’) 

%%% 

fn  =  fn  +  1 ;  figure(fn),  orient  tall 
plot(SNR_est(:,:,fn),  BER(:,:,fn),  'b.') 
hold  on 

plot(mean(SNR_est(:,:,fn),2),  mean(BER(:,:,fn),2),  'r-') 
hold  off 

title(['SNR  =  ’,num2str(SNR),  ’dB,  f_c  =  ’,num2str(fc/1000),'kHz,  R 

',num2str(R/l 000), 'kbps,  Range  =  4  km,  Rough  =  4,  Water  Depth  =  100  m’]) 

xlabel('SNR  (dB)’),  ylabel(’BER’) 
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%%% 

fc  =  6e3 

R  =  (5+0/3)*  le3 

%%% 

fn  =  fn  +  1 ;  figure(fn),  orient  tall 
plot(SNR_est(:,:,fn),  BER(:,:,fn),  ’b.') 
hold  on 

plot(mean(SNR_est(:,:,fn),2),  mean(BER(:,:,fn),2),  ’r-’) 
hold  off 

title(['SNR  =  ’,num2str(SNR),  ’dB,  f_c  =  ',num2str(fc/1000),'kHz,  R 
',num2str(R/l  000), 'kbps,  Range  =  2  km,  Rough  =  0,  Water  Depth  =  340  m']) 

xlabel('SNR  (dB)’),  ylabel(’BER’) 

%%% 

fn  =  fn  +  1 ;  figure(fn),  orient  tall 
plot(SNR_est(:,:,fn),  BER(:,:,fn),  'b.') 
hold  on 

plot(mean(SNR_est(:,:,fn),2),  mean(BER(:,:,fn),2),  ’r-’) 
hold  off 

title(['SNR  =  ’,num2str(SNR),  'dB,  f_c  =  ',num2str(fc/1000),'kHz,  R 

',num2str(R/l 000), 'kbps,  Range  =  2  km,  Rough  =  2,  Water  Depth  =  340  m']) 

xlabel('SNR  (dB)’),  ylabel(’BER’) 

%%% 

fn  =  fn  +  1 ;  figure(fn),  orient  tall 
plot(SNR_est(:,:,fn),  BER(:,:,fn),  'b.') 
hold  on 

plot(mean(SNR_est(:,:,fn),2),  mean(BER(:,:,fn),2),  'r-') 
hold  off 

title(['SNR  =  ’,num2str(SNR),  'dB,  f_c  =  ’,num2str(fc/1000),'kHz,  R 

’,num2str(R/l 000), ’kbps,  Range  =  4  km,  Rough  =  0,  Water  Depth  =  340  m’]) 

xlabel(’SNR  (dB)’),  ylabel(’BER’) 

%%% 

fn  =  fn  +  1 ;  figure(fn),  orient  tall 
plot(SNR_est(:,:,fn),  BER(:,:,fn),  ’b.’) 
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hold  on 

plot(mean(SNR_est(:,:,fn),2),  mean(BER(:,:,fn),2),  'r-') 
hold  off 

title(['SNR  =  ’,num2str(SNR),  ’dB,  f_c  =  ',num2str(fc/1000),'kHz,  R  = 
’,num2str(R/l  000), 'kbps,  Range  =  4  km,  Rough  =  2,  Water  Depth  =  340  m']) 

xlabel('SNR  (dB)’),  ylabel(’BER’) 

%%% 

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 


%64ary  QAM  PERFORMANCE 
SNR  =  27 
fc  =  8e3 

R  =  (l  0+0/3)*  le3 
%%% 

fn  =  fin  +  1 ;  figure(fn),  orient  tall 
plot(SNR_est(:,:,fn),  BER(:,:,fn),  'b.') 
hold  on 

plot(mean(SNR_est(:,:,fn),2),  mean(BER(:,:,fn),2),  'r-') 
hold  off 

title(['SNR  =  ’,num2str(SNR),  'dB,  f_c  =  ',num2str(fc/1000),'kHz,  R  = 

’,num2str(R/l  000), 'kbps,  Range  =  2  km,  Rough  =  0,  Water  Depth  =  100  m']) 

xlabel('SNR  (dB)’),  ylabel(’BER’) 

%%% 

fn  =  fn  +  1 ;  figure(fn),  orient  tall 
plot(SNR_est(:,:,fn),  BER(:,:,fn),  'b.') 
hold  on 

plot(mean(SNR_est(:,:,fn),2),  mean(BER(:,:,fn),2),  'r-') 
hold  off 

title(['SNR  =  ’,num2str(SNR),  'dB,  f_c  =  ',num2str(fc/1000),'kHz,  R  = 

',num2str(R/l  000), 'kbps,  Range  =  4  km,  Rough  =  0,  Water  Depth  =  100  m']) 

xlabel('SNR  (dB)’),  ylabel(’BER’) 

%%% 

fc  =  6e3 
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R=  (5+0/3)*  le3 

%%% 

fn  =  fn  +  1 ;  figure(fn),  orient  tall 
plot(SNR_est(:,:,fn),  BER(:,:,fn),  'b.’) 
hold  on 

plot(mean(SNR_est(:,:,fn),2),  mean(BER(:,:,fn),2),  'r-') 
hold  off 

title(['SNR  =  ’,num2str(SNR),  ’dB,  f_c  =  ’,num2str(fc/1000),'kHz,  R 
’,num2str(R/l 000), ’kbps,  Range  =  2  km,  Rough  =  0,  Water  Depth  =  100  m’]) 

xlabel(’SNR  (dB)'),  ylabel(’BER’) 

%%% 

fn  =  fn  +  1 ;  figure(fn),  orient  tall 
plot(SNR_est(:,:,fn),  BER(:,:,fn),  'b.') 
hold  on 

plot(mean(SNR_est(:,:,fn),2),  mean(BER(:,:,fn),2),  ’r-’) 
hold  off 

title(['SNR  =  ’,num2str(SNR),  ’dB,  f_c  =  ’,num2str(fc/1000),'kHz,  R 

’,num2str(R/l  000), 'kbps,  Range  =  2  km,  Rough  =  2,  Water  Depth  =  100  m']) 

xlabel('SNR  (dB)’),  ylabel(’BER’) 

%%% 

fn  =  fn  +  1 ;  figure(fn),  orient  tall 
plot(SNR_est(:,:,fn),  BER(:,:,fn),  ’b.’) 
hold  on 

plot(mean(SNR_est(:,:,fn),2),  mean(BER(:,:,fn),2),  ’r-’) 
hold  off 

title(['SNR  =  ’,num2str(SNR),  'dB,  f_c  =  ',num2str(fc/1000),'kHz,  R 

’,num2str(R/l 000), 'kbps,  Range  =  2  km,  Rough  =  4,  Water  Depth  =  100  m']) 

xlabel('SNR  (dB)’),  ylabel(’BER’) 

%%% 

fn  =  fn  +  1 ;  figure(fn),  orient  tall 
plot(SNR_est(:,:,fn),  BER(:,:,fn),  'b.') 
hold  on 

plot(mean(SNR_est(:,:,fn),2),  mean(BER(:,:,fn),2),  'r-') 
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hold  off 


title(['SNR  =  ’,num2str(SNR),  ’dB,  f_c  =  ’,num2str(fc/1000),'kHz,  R 
',num2str(R/l  000), 'kbps,  Range  =  4  km,  Rough  =  0,  Water  Depth  =  100  m']) 

xlabel('SNR  (dB)’),  ylabel(’BER’) 

%%% 

fn  =  fn  +  1 ;  figure(fn),  orient  tall 
plot(SNR_est(:,:,fn),  BER(:,:,fn),  'b.') 
hold  on 

plot(mean(SNR_est(:,:,fn),2),  mean(BER(:,:,fn),2),  ’r-’) 
hold  off 

title(['SNR  =  ’,num2str(SNR),  'dB,  f_c  =  ',num2str(fc/1000),'kHz,  R 

’,num2str(R/l 000), 'kbps,  Range  =  4  km,  Rough  =  2,  Water  Depth  =  100  m']) 

xlabel('SNR  (dB)’),  ylabel(’BER’) 

%%% 

fn  =  fn  +  1 ;  figure(fn),  orient  tall 
plot(SNR_est(:,:,fn),  BER(:,:,fn),  'b.') 
hold  on 

plot(mean(SNR_est(:,:,fn),2),  mean(BER(:,:,fn),2),  'r-') 
hold  off 

title(['SNR  =  ’,num2str(SNR),  'dB,  f_c  =  ’,num2str(fc/1000),’kHz,  R 

’,num2str(R/l  000), 'kbps,  Range  =  4  km,  Rough  =  4,  Water  Depth  =  100  m']) 

xlabel('SNR  (dB)’),  ylabel(’BER’) 

%%% 

fc  =  6e3 

R  =  (5+0/3)*  1  e3 

%%% 

fn  =  fn  +  1 ;  figure(fn),  orient  tall 
plot(SNR_est(:,:,fn),  BER(:,:,fn),  'b.') 
hold  on 

plot(mean(SNR_est(:,:,fn),2),  mean(BER(:,:,fn),2),  'r-') 
hold  off 

title(['SNR  =  ’,num2str(SNR),  'dB,  f_c  =  ’,num2str(fc/1000),'kHz,  R 
',num2str(R/l 000), 'kbps,  Range  =  2  km,  Rough  =  0,  Water  Depth  =  340  m']) 

xlabel('SNR  (dB)’),  ylabel(’BER’) 
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%%% 

fn  =  fn  +  1 ;  figure(fn),  orient  tall 
plot(SNR_est(:,:,fn),  BER(:,:,fn),  ’b.’) 
hold  on 

plot(mean(SNR_est(:,:,fn),2),  mean(BER(:,:,fn),2),  'r-') 
hold  off 

title(['SNR  =  ’,num2str(SNR),  ’dB,  f_c  =  ’,num2str(fc/1000),'kHz,  R 
’,num2str(R/l  000), 'kbps,  Range  =  2  km,  Rough  =  2,  Water  Depth  =  340  m']) 

xlabel('SNR  (dB)’),  ylabel(’BER’) 

%%% 

fn  =  fn  +  1 ;  figure(fn),  orient  tall 
plot(SNR_est(:,:,fn),  BER(:,:,fn),  'b.') 
hold  on 

plot(mean(SNR_est(:,:,fn),2),  mean(BER(:,:,fn),2),  'r-') 
hold  off 

title(['SNR  =  ’,num2str(SNR),  'dB,  f_c  =  ',num2str(fc/1000),'kHz,  R 

',num2str(R/l  000), 'kbps,  Range  =  4  km,  Rough  =  0,  Water  Depth  =  340  m']) 

xlabel('SNR  (dB)’),  ylabel(’BER’) 

%%% 

fn  =  fn  +  1 ;  figure(fn),  orient  tall 
plot(SNR_est(:,:,fn),  BER(:,:,fn),  'b.') 
hold  on 

plot(mean(SNR_est(:,:,fn),2),  mean(BER(:,:,fn),2),  'r-') 
hold  off 

title(['SNR  =  ’,num2str(SNR),  'dB,  f_c  =  ’,num2str(fc/1000),'kHz,  R 

',num2str(R/l 000), 'kbps,  Range  =  4  km,  Rough  =  2,  Water  Depth  =  340  m']) 

xlabel('SNR  (dB)’),  ylabel(’BER’) 
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